Abstract
Background
Molecular evolution is usually described assuming a neutral or weakly nonneutral substitution model. Recently, new data have become available on evolution of sequence regions under a selective pressure, e.g. transcription factor binding sites. To reconstruct the evolutionary history of such sequences, one needs evolutionary models that take into account a substantial constant selective pressure.
Results
We present a simple evolutionary model with a single preferred (consensus) nucleotide and the neutral substitution model adopted for all other nucleotides. This evolutionary model has a rate matrix in which all substitutions that do not involve the consensus nucleotide occur with the same rate. The model has two time scales for achieving a stationary distribution; in the general case only one of the two rate parameters can be evaluated from the stationary distribution. In the middletime zone, a counterintuitive behavior was observed for some parameter values, with a probability of conservation for a nonconsensus nucleotide greater than that for the consensus nucleotide. Such an effect can be observed only in the case of weak preference for the consensus nucleotide, when the probability to observe the consensus nucleotide in the stationary distribution is less than 1/2. If the substitution rate is represented as a product of mutation and fixation, only the fixation can be calculated from the stationary distribution. The exhibited conservation of nonconsensus nucleotides does not take place if the elements of mutation matrix are identical, and can be related to the reduced mutation rate between the nonconsensus nucleotides. This bias can have no effect on the stationary distribution of nucleotide frequencies calculated over the ensemble of multiple alignments, e.g. transcription factor binding sites upstream of different sets of coregulated orthologous genes.
Conclusion
The derived model can be used as a null model when analyzing the evolution of orthologous transcription factor binding sites. In particular, our findings show that a nucleotide preferred at some position of a multiple alignment of binding sites for some transcription factor in the same genome is not necessarily the most conserved nucleotide in an alignment of orthologous sites from different species. However, this effect can take place only in the case of a mutation matrix whose elements are not identical.
Background
The controlled expression of genes is the main mechanism responsible for the life cycle and biodiversity [1]. Transcription, which is a crucial process defining the level of gene expression, is regulated by interaction of transcription factor proteins (TFs) with transcription factor binding sites (TFBSs) in a DNA molecule [2]. Thus, adaptable interaction between TFs and TFBSs is one of the main driving forces of biological evolution [3].
Both TFs and TFBSs are subject to mutations and selection affecting their interaction [4]. In this study we focus on mutations in TFBSs. New experimental [57] and computational [8] methods of TFBS identification produce an increasing amount of data about TFBS sequences, which creates a possibility to study evolutionary events in these regions.
Modelling evolution of regulatory sequences can be useful for understanding both the general mechanisms of gene expression control and the regulatory history of particular genes.
The evolution of regulatory regions has a complex pattern [3,911] and is still unclear in many aspects. DNA segments can gain and lose the TFBS function [12], and that can bring new genes under regulation by a particular TF [1315], or divert regulation from other genes [16,17]. One particular type of events is emergence of new or changed sites following displacement of a transcription factor by horizontal transfer [18]. Sometimes this leads to considerable changes in the regulon content [19,20] or even partial or complete rewiring of regulatory cascades [2126], reviewed in [27].
Evolution of functional TFBS sequences is strongly nonneutral [11,28] and under a positive selection [29], which makes it difficult to calculate the rate of TFBS evolution. This rate varies between TFBS positions [30,31]. Moreover, coevolution of TFBS and TF [16,20] can make the selective pressure vary in different lineages. Indeed, although in some cases the DNA motif bound by orthologous factors may be conserved at surprisingly large evolutionary distances [14,3234], in other cases not only the motifs themselves may be different [35,36], but even the symmetry of the motif (e.g. palindrome or direct repeat) may change [37,38].
The existing evolutionary models, which were successful in reconstruction of phylogenetic relations, can be applied to evolution of regulatory sequences only with a caution. Such models are historically related to the Jukes–Cantor [39] and Kimura [40] models of molecular evolution. Existing modifications of these models take into account various global characteristics like transition/transversion rate or local GC composition [4144]. They are not applicable to the case of strong selection for a specific nucleotide at a particular position.
On the other hand, models developed specifically for the evolution of TFBS are needed to reconstruct the evolutionary origin of a particular TFBS and to evaluate the positionspecific mutation rate and selective pressure.
Because of the positionspecific variation in the rate of TFBS evolution [30], the rate matrix must also be positionspecific. The data produced by mass experiments on TFBS identification or comparative genomic studies produce tens of TFBS for each TF (more exactly, a group of orthologous TFs). That might be sufficient to evaluate the evolutionary rate at each TFBS position.
Here we consider the simplest model of positionspecific evolution with one preferred (consensus) nucleotide and three other (minor) nucleotides, the latter considered in a symmetric setting, without any selection or rate preferences [45]. Such a model can be deduced from physical requirements of the TF/TFBS interaction [46] and can explain the observed TFBS fuzziness.
We build a rate matrix, which enhances the model of [45]. We calculate the substitution probability for each finite time and show that the nucleotide conservation in phylogenetic lineages can be nontrivial for some parameter values. Particularly, a nonconsensus nucleotide may appear more conserved than the consensus nucleotide, although the latter has a selective preference. This happens when the rate of mutations between nonconsensus nucleotides is lower than the rate of mutation into the consensus, or there is selection against any mutation in nonconsensus nucleotides.
Results and Discussion
Model
We start with definitions. We consider an alignment of several sequences; all positions in this alignment are assumed to be independent, and thus may be modelled independently. Consensus nucleotide (or simply consensus) is the most frequent nucleotide in an alignment column (position). Other nucleotides are called nonconsensus. The frequency of the consensus nucleotide N_{c }is a fraction of the number of consensus nucleotides in a particular position. Obviously, 1/4 <N_{c }≤ 1. The consensus is called weak, if 1/4 <N_{c }< 1/2; the consensus is strong, if 1/2 ≤ N_{c }< 1.
Consider the model of nucleotide substitutions given by a Markov process X(t) with four states {g_{1}, g_{2}, g_{3}, g_{4}}. Without loss of generality, assume that the state g_{1 }is the consensus state and the states g_{2}, g_{3}, g_{4 }are equiprobable nonconsensus states. Suppose that the transition rate matrix A = (q_{ij}) is given by
where q_{ij }is the transition rate from the state g_{i }to g_{j }and α, β, and δ are positive unknown parameters.
Transitional probabilities
Let P_{ij}(t) be the transitional probability from the state g_{i }to the state g_{j }for the time t. From the theory of Markov chains we have for h → 0
For brevity we denote by 'c' the subscript '1' corresponding to the consensus state g_{1 }and by 'n' and 'm' the subscripts 2, 3, and 4 corresponding to the nonconsensus states g_{2}, g_{3}, g_{4}. Thus we have five transitional probabilities P_{cc}, P_{cn}, P_{nc}, P_{nn}, and P_{nm}, where n and m stand for two distinct nonconsensus states. The formulas for P_{ij }are derived from simple calculation:
For simplicity, introduce a new timescale, u = tβ, and denote λ = α/β, μ = δ/β. Then
and
Note that the parameter μ is related to the conservation within the set of nonconsensus states {g_{2}, g_{3}, g_{4}}. Simply, λ is a rate of transition from the consensus nucleotide and μ is a rate of transition between nonconsensus states up to the timescale parameter β.
The stationary distribution π = (π_{c}, π_{n}, π_{n}, π_{n}) of the process X(t) is given by
According to our definition of the consensus as the most frequent nucleotide, from π_{c }= 1/(3λ + 1) > 1/4 it follows that 0 <λ < 1. The parameter λ is responsible for transitions from the consensus state to the nonconsensus ones. At the same time, the transition from a nonconsensus state to the consensus state occurs with the rate 1 (up to the timescale parameter β). Intuitively, in the model with one selected consensus state the probability of transition to a nonconsensus state should be smaller than the probability of transition to the consensus state; thus, λ should be less than 1. Indeed, if λ ≥ 1, then π_{c }≤ 1/4, π_{n }≥ 1/4 and we have three equiprobable consensus states.
If the consensus is strong, then π_{c }≥ 1/2, and, consequently, 0 <λ ≤ 1/3. Note that we can not impose any condition on μ other than μ > 0.
Conservation and transition of nucleotides
The goal of this section is to compare the transitional probabilities between different states in our model. We will see that the relations we get have clear biological interpretation.
Simple cases
Let a unique consensus exist in our model, i.e. 0 <λ < 1. It is not difficult to obtain the following relations between the transitional probabilities, which hold for any λ ∈ (0, 1), u > 0.
1. P_{cc}(u) > P_{cn}(u);
2. P_{cc}(u) > P_{nc}(u);
3. P_{cc}(u) > P_{nm}(u) for any μ > 0;
4. P_{nn}(u) > P_{cn}(u) for any positive λ and μ;
5. P_{nc}(u) > P_{cn}(u);
6. P_{nn}(u) > P_{nm}(u) for any μ > 0.
The interpretation of these results is straightforward. Figures 1, 2, 3, 4 show the graphs of the transitional probabilities for different values of λ and μ. The relations 1–6 between the transitional probabilities are clearly shown in the figures. The first three inequalities show that the probability of conservation of the consensus state is always higher than the probability of transition between the states of different type. Inequalities 1, 4, and 5 imply that the probability of transition from the consensus state to a nonconsensus one is always less than the probability of conservation of a nucleotide or the probability of transition from a nonconsensus state to the consensus one. Inequalities 3 and 6 show that the probability of transition between two different nonconsensus states is always less than the probability of state conservation. All these results correlate well with our intuitive ideas about an evolutionary model with selective pressure. Note that relations 1–6 hold for any positive μ, λ ∈ (0, 1) and for any time period.
Figure 1. Transitional probabilities P_{cn}, P_{cc}, P_{nn}, P_{nc}, P_{nm }for λ = 1/6, μ = 1/2. This figure describes the case P_{cc}(u) > P_{nn}(u) for all u > 0. The consensus is strong, since π_{c }= 2/3. Since the rate of transition between nonconsensus states μ is greater than the rate of transition from the consensus state λ, we have P_{cn}(u) <P_{nm}(u). The relation between the probability of conservation of a nonconsensus nucleotide and the probability of transition from the nonconsensus state to the consensus one is shown, P_{nn}(u) > P_{nc}(u) for u ∈ (0, u_{0}) and, vice versa, P_{nn}(u) <P_{nc}(u) for u > u_{0}. This figure exemplifies all simple cases that hold for any λ ∈ (0, 1), μ > 0, u > 0 (see the corresponding paragraph in the main text).
Figure 2. Transitional probabilities P_{cn}, P_{cc}, P_{nn}, P_{nc}, P_{nm }for λ = 2/3, μ = 3/4. Here we have P_{nn}(u) > P_{nc}(u) for u ∈ (0, u_{0}) and, vice versa, P_{nn}(u) <P_{nc}(u) for u > u_{0}. In this figure P_{cc}(u) > P_{nn}(u) for all u > 0, as 3λ ≤ 2μ + 1. The consensus is weak, since π_{c }= 1/3. Since μ > λ, we have P_{cn}(u) <P_{nm}(u). This figure shows all simple cases that hold for any λ ∈ (0, 1), μ > 0, u > 0.
Figure 3. Transitional probabilities P_{cn}, P_{cc}, P_{nn}, P_{nc}, P_{nm }for λ = 2/3, μ = 1/6. In this case the probability of conservation of the consensus nucleotide is less than the probability of conservation of a nonconsensus nucleotide, P_{cc}(u) <P_{nn}(u) for u ∈ (0, u**), since 3λ > 2μ + 1 and λ > μ. The stationary distribution of the consensus state is π_{c }= 1/3. Thus, the consensus is weak. Since μ <λ, we have P_{cn}(u) > P_{nm}(u) for all u > 0.
Figure 4. Transitional probabilities P_{cn}, P_{cc}, P_{nn}, P_{nc}, P_{nm }for λ = 1/6, μ = 2. Here we have P_{cc}(u) > P_{nn}(u) for all u > 0, since 3λ ≤ 2μ +1. The consensus is strong and the stationary distribution of the consensus state is π_{c }= 2/3. Since μ > λ, we have P_{cn}(u) <P_{nm}(u). The probability of substitutions between nonconsensus nucleotides P_{nm}(u) is not monotone with a maximum point at u_{2 }= 2 log 7/11 ≈ 0.354. Since μ > 1, transitions between nonconsensus nucleotides occur more frequently than transitions from a nonconsensus to the consensus nucleotide, P_{nc}(u) <P_{nm}(u) on the time interval u ∈ (0, u_{1}), whereas P_{nc}(u) > P_{nm}(u) for u > u_{1}. This figure also shows the relation between P_{nn}(u) and P_{nc}(u).
Interesting cases
The most interesting case is the conservation of the consensus or a nonconsensus nucleotide. Recall that P_{cc}(h) = 1  3 λβ h + o(1), P_{nn}(h) = 1  (2μ + 1)β h + o(1) for h → 0. Thus, the relation between the probabilities of conservation of the consensus state and conservation of a nonconsensus state during the time u depends on the relation between 3λ and 2μ + 1. We obtained the following result.
1. If 3λ ≤ 2μ + 1, then
2. If 3λ > 2μ + 1, λ > μ, then
Here u** = u**(λ, μ) > 0 is a time moment depending on λ and μ. More precisely, u** is a nonzero solution of the equation P_{cc}(u) = P_{nn}(u). Let F (u) = P_{cc}(u)  P_{nn}(u). It can be shown that for λ > (2μ + 1)/3 and λ > μ
Indeed,
and F'(u) = 0 if and only if
2(3μ + 1) e^{(3μ+1)u }= (9λ  1)e^{(3λ+1)u}.
The solution of the latter equation is u*. If 3λ > 2μ + 1 and λ > μ, then u* > 0 and F increases for u > u* and decreases for u <u*. At the same time, F(0) = 0. Thus, in this case F(u) < 0 for u ∈ (0, u**) and F(u) > 0 for u > u**. This implies the second statement.
Further, if 3λ ≤ 2μ + 1 and λ > μ, then u* < 0 and F'(u) > 0 for all u > 0. Therefore, F(u) increases for u > 0 and F(u) > 0, since F(0) = 0. Next, for λ = μ we obtain that
for any u > 0. Thus, F(u) > 0 if 3λ ≤ 2μ + 1, and the first inequality follows.
The second relation partly describes the case of the weak consensus (Fig. 3). Since 3λ + 1 > 2(μ + 1), the stationary distribution of the consensus state is estimated from above as
At the same time, the first relation holds both in the case of strong consensus π_{c}≥ 1/2 (Fig. 1) and in the case of weak consensus 1/4 <π_{c }< 1/2 (Fig. 2).
The second interesting result is a relation between the probability of conservation of the nonconsensus nucleotide and the probability of transition from a nonconsensus state to the consensus state. We have
where u_{0 }≡ u_{0}(λ, μ) is the solution of equation P_{nn}(u) = P_{nc}(u). It can be shown that u_{0 }is always positive. This relation shows that the probability of conservation of a nonconsensus nucleotide is less than the probability of transition to the consensus nucleotide from a nonconsensus one for sufficiently long time u > u_{0}. However, on the interval (0, u_{0}) the opposite inequality holds (see Fig. 1, 2).
Less interesting cases
It remains to study the relations between P_{nm }and P_{cn}, P_{nc}. From the practical point of view, these cases are not very interesting, since the events of transition between different nonconsensus states can be hardly observed on practice. However, we consider them for the sake of completeness.
The following relations imply that the transitional probability from the consensus to a nonconsensus state could be less or greater than the transitional probability between two different nonconsensus states depending on the relation between μ and λ (see Fig. 2 and 3):
Indeed, if μ > λ, then δ > α and the transition rate d between the nonconsensus states is greater than the transition rate from the consensus state to nonconsensus. Clearly, in this case P_{nm }has to be greater than P_{cn }(see Fig. 2 and Fig. 3).
The next case concerns the relation P_{nm }and P_{nc }for λ ∈ (0, 1) (see Fig. 1 and 4):
where u_{1 }= u_{1}(λ, μ) is a solution of the equation P_{nc}(u) = P_{nm}(u).
An interesting fact is that the probability P_{nm}(u) is not monotone in u for μ > λ (see Fig. 4). If μ > λ, then P_{nm }increases on (0, u_{2}) and decreases for u > u_{2}, where
Finally, consider the case μ = λ. As it has been shown above, P_{cc}(u) > P_{nn}(u) for u > 0 (see Fig. 5). It is easy to see that in this case P_{nm }≡ P_{cn}. The function P_{nm}(u) is monotonically increasing for all u. Note that in the degenerate case λ = μ = 1 all states are equiprobable, the stationary distribution of the process π = (1/4, 1/4, 1/4, 1/4) and P_{nn }= P_{cc}, P_{nc }= P_{cn }= P_{nm }(see Fig. 6).
Figure 5. Transitional probabilities for λ = μ = 1/2. This case is similar to the one of Fig. 1, since here 3λ ≤ 2μ + 1 and, consequently, the probability of consensus conservation is always greater than the probability of conservation of a nonconsensus nucleotide, P_{cc}(u) > P_{nn}(u) for all u > 0. As λ = μ, we have P_{nm }= P_{cn}.
Figure 6. Degenerate case λ = μ = 1. In this case all states are equiprobable, the stationary distribution of the process π = (1/4, 1/4, 1/4, 1/4) and P_{nn }≡ P_{cc}, P_{nc }≡ P_{cn }≡ P_{nm}.
Generalization
The obtained results can be generalized on the following model. Consider a Markov chain (t) with M + N states {g_{1}, ..., g_{M}, g_{M+1}, ..., g_{M+N }}, where the states g_{1}, ..., g_{M }are consensus states, and the states g_{M+1}, ..., g_{M+N }are nonconsensus ones. Define the transition rate matrix for this process by
where α, β, γ, δ, are positive unknown parameters.
In this case we have two groups of M consensus and N nonconsensus states. We use the same notation for subscripts of transitional probabilities . If c and d denote different consensus states, n and m stand for different nonconsensus states as before, we have the following transitional probabilities:
Thus, this process has the stationary distribution = (π_{c}, ..., π_{c}, π_{n}, ..., π_{n}) with
Clearly, β > γ, since the first M states are the consensus ones. Next, compare the probabilities of conservation of the consensus and nonconsensus states P_{cc }and P_{nn}. In a similar way, we obtain that there exists t* > 0 such that P_{cc}(t) <P_{nn}(t) for t ∈ (0, t*). Analyzing the probabilities P_{cc }and P_{nn }we can show that this is possible only for δ (N  1)  α (M  1) + β M  γ N < 0. Then the frequency of the consensus state (stationary distribution of consensus) is estimated from above as
If M = 1, N = 3, then
This condition coincides with the condition on weak consensus obtained for the case of four states with one consensus (M = 1, N = 3) considered above.
Comparison with the Molecular Evolution Theory
In the framework of the molecular evolution theory, the element a_{ij }of the transition rate matrix is considered to be proportional to the product of the mutation rate p_{ij }and the probability of fixation of a mutation f_{ij}, a_{ij }= kp_{ij}f_{ij}, where k is an arbitrary scaling constant [47]. As in [45] we start from the simplest Jukes–Cantor (1  p)scheme, to which we introduce selection. Thus, we ignore the difference between transitions and transversions in p_{ij}.
TFBS regulating different genes in the same genome most likely evolve independently and thus the nucleotide composition π_{i }at the respective positions of different TFBS occurrences approximates the equilibrium frequencies. With these equilibrium frequencies at hand it is possible to relate p_{ij }and f_{ij }with the equation (see [47])
In our case, for the substitution rate between three nonconsensus positions we obtain π_{i }= π_{j }= π_{n }and p_{ij }= p_{ji }= p_{nm}, which yields f_{nm}∝ 1 by the l'Hôpital rule as in [47]. Thus, δ = r_{nm }= k_{pnm}.
For the substitutions between nonconsensus and consensus positions, r_{nc}, both the selection preferences and mutation asymmetry come into consideration. In this case the "asymmetry constant" λ is crucial, which satisfies the inequality r_{cn}/r_{nc }= p_{n}/p_{c }= λ < 1. The following expression is valid:
This fixation rate is linked with the Kimura selection constant [48] s by the relation f_{cn }= (1  e^{2s})/(1  e^{2Ns}), where N is the population size.
If the mutation rate is symmetric, p_{nc }= p_{cn}, then f_{cn }∝ λ log(1/λ)/(1  λ). Conversely, for the nonconsensus to consensus substitution f_{nc}∝ log(1/λ)/(1  λ) = f_{cn}/λ, the greater flux from a nonconsensus state to the consensus maintains a greater consensus frequency. Note that in alignments of sites in a single genome we can observe only the equilibrium constants π_{n}, π_{c }(which actually are rather rough approximations), thus the assumption p_{nc }= p_{cn }may be too strong, and the above general formula for f_{cn }might be more relevant.
The coefficients in the matrix A for the symmetric mutation rates are given by
If the background mutation rate is identical for all consensus and nonconsensus nucleotides, we obtain p_{cn }= p_{nc }= p_{nm }= p and p may be merged with the constant k. In this most simple case we obtain
and, consequently, β > δ > α . This is the simplest generalization of the Jukes–Cantor model for the case with introduced selection.
It should be noted that in this case μ = δ/β = (1  λ)/log(1/λ), which implies μ > λ. Thus, 3λ ≤ 2μ + 1 and we are in Case 1 of "Interesting Cases" for which P_{cc}(u) > P_{nn}(u) for u > 0.
Conservation of nonconsensus nucleotides at the reduced mutation rate
Previously [8] we have observed that nonconsensus nucleotides may be highly conserved in the alignments of orthologous TFBS in bacterial genomes. On the other hand, as shown above, the nonconsensus nucleotide in the alignment of orthologous sites from different species cannot be more conserved than the consensus nucleotide if we adopt the Kimura model with the identical mutation rate for all pairs.
One way to explain the observation made in [8] is to drop the equivalence of nonconsensus nucleotides and assume that different nonconsensus nucleotides are under different selection in the sites regulating different rows of orthologous genes. Interestingly, it is possible to observe such conservation pattern in the model with an identical probabilities of mutation and fixation for all nonconsensus nucleotides in all orthologous site rows, with a single preferred nucleotide, the consensus, the same in all sites. The only necessary relaxation of the model is to drop the condition p_{cn }= p_{nc }= p_{nm }= p for the condition p_{nm }<p_{nc }= p_{cn}. Doing this it is possible to satisfy the inequalities determining Case 2 of "Interesting Cases": μ <λ, 3λ > 2μ + 1.
The condition p_{nm }<p_{nc }= p_{cn }means that the rate of direct mutations from one nonconsensus nucleotide to another one is lower than the mutation rate in pairs involving the consensus nucleotide. A possible example of such specific reduction of the mutation rate comes from correlations of nucleotides occupying different positions of the same binding site. Assume that two positions within the site are not independent and must be occupied by correlated nucleotides. These may be, e.g., adjacent positions in a DNA site or basepaired positions in an RNA structure. Assume also that if any of the two positions is occupied by the consensus nucleotide, the correlating nucleotide may be arbitrary. Conversely, if one position is occupied by a nonconsensus nucleotide, the other position should be occupied with some specific nucleotide, e.g. the complementary one in the case of an RNA structure.
In this case, the preferred pathway from a nonconsensus nucleotide to another nonconsensus nucleotide would become not via a direct mutation, but via an intermediate mutation into the consensus nucleotide. For example, if "C" is both cytosine and consensus, and for some case this position is correlated with another one, so that only AA, TT, and GG pairs involving the nonconsensus nucleotides at the first position are allowed, then only mutation A>C is valid, whereas mutations A>T and A>G are forbidden, and may occur via mutation into C and then the compensating mutation in the second position of the site. This simple model to some extent agrees with recent studies demonstrating that proteinDNA interactions are rather complex and probably may not be described by a simple positionindependent model such as a positional weight matrix [49].
At the same time, this effect is not in contradiction with the uniform distribution of nonconsensus nucleotides obtained from alignments of multiple TFBS regulating different genes in the same genome. It also allows for high conservation of some nonconsensus nucleotides in the alignment of orthologous transcription factor binding sites from different species. Indeed, if the contextdependent pattern of conservation is specific to a particular position in a particular set of orthologous TFBSs, exactly this type of behavior may be expected.
Conclusion
The evolutionary model derived here can be generalized to the case of M consensus and N nonconsensus nucleotides that are equiprobable, respectively. Thus, for M = 2, N = 2 and α = β, γ = δ. we get the case of neutral evolution with different rates of transitions and transversions [40]. If M = 1, N = 3, then we have the evolution model under constant selective pressure considered in this paper.
One somewhat nonobvious feature of this model is the existence of a combination of the rate of transition between nonconsensus states μ and the rate of transition from the consensus state λ, for which the conservation of the consensus nucleotide in a sequence alignment can be lower than the conservation of a neutral (nonconsensus) nucleotide. However, this can be observed only for the case of a weak consensus 1/4 <N_{c }< 1/2 and for a relatively short time interval u ∈ (0, u**), and a mutation matrix with different elements, e.g. contextdependent.
Models of this type can be applied not only for the analysis of regulatory sites, but in other situations, e.g. for the analysis of functional sites in proteins [50] or analysis of evolution in the case of nucleotide biases [4144,51,52].
In future work, we intend to estimate the parameters of the model from the data based on the method of maximum likelihood trees. Consider a Markov process X(t) describing the evolution at some fixed position. We assume that we observe only the endpoints of k different paths of X(t) depending on the evolution branch. In other words, our data is a set of k nucleotides at a fixed position in k genomes. The parameter λ can be easily estimated by the 'naive' estimator = (1/N_{c } 1)/3, where N_{c }is the frequency of the consensus nucleotide at the fixed position within k genomes. This estimator is obtained from the formula for the stationary distribution of the consensus state π_{c }= 1/(3λ + 1). However, it is a good approximation of the true value of λ. The preliminary analysis of numerically simulated data shows that the performance of the maximum likelihood estimator is also rather good. On the other hand, it is not clear whether it is possible to construct an explicit estimator for the rate of transition between nonconsensus states μ from the data, although μ participates in the expression for transitional probabilities.
The inspiration for the constructed model was the example of a set of orthologous transcription factors interacting with the cognate regulatory regions. However, it appears that evolution of sequences under a low selection pressure is a more widespread phenomenon. Indeed, a recent study [53] demonstrated that the force causing conservation of some noncoding genome regions of human, mouse and chimpanzee can be explained by a rather small selective pressure at the genomic level. A similar problem appears in the context of the CG composition of genomic regions [54]. Again, in this case the selective pressure appears to be low, although unlike the previous examples, here two nucleotides become selected for rather than one consensus nucleotide. The appropriate model in this case includes differences in the transition and the transversion rates, which makes the model more complicated, and probably would result in more complex time behavior of the substitution probabilities.
Anyhow, the emerging huge amount of data on orthologous noncoding regions, which has became available recently, brings forward a problem of modelling evolution with a selection pressure at finite times.
Methods
The numerical simulations and figures were produced using Matlab.
Authors' contributions
MSG and VJM conceived the study. FNE and EAK developed the model. FNE performed numerical simulations. FNE, VJM and MSG wrote the paper. All authors have read and approved the final version.
Acknowledgements
This study was partially supported by the Russian Federation Agency in Science and Innovation State Contract # 02.531.11.9003 (VJM, EAK), by the Russian Foundation for Basic Research, projects no. 070401584 (VJM), 070401623 (EAK) and by grants from the Howard Hughes Medical Institute (55005610, MSG), INTAS (0510000088028, VJM and MSG), and the Russian Academy of Sciences (program "Molecular and Cellular Biology", MSG). FNE was supported by the Russian Science Support Fund.
References

Levine M, Tjian R: Transcription regulation and animal diversity.
Nature 2003, 424(6945):147151. PubMed Abstract  Publisher Full Text

Pabo CO, Sauer RT: Transcription factors: structural families and principles of DNA recognition.
Annu Rev Biochem 1992, 61:10531095. PubMed Abstract  Publisher Full Text

Wray GA, Hahn MW, Abouheif E, Balhoff JP, Pizer M, Rockman MV, Romano LA: The evolution of transcriptional regulation in eukaryotes.
Mol Biol Evol 2003, 20(9):13771419. PubMed Abstract  Publisher Full Text

Mirny LA, Gelfand MS: Structural analysis of conserved base pairs in proteinDNA complexes.
Nucleic Acids Res 2002, 30(7):17041711. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Wells J, Farnham PJ: Characterizing transcription factor binding sites using formaldehyde crosslinking and immunoprecipitation.
Methods 2002, 26:4856. PubMed Abstract  Publisher Full Text

Weinmann AS, Farnham PJ: Identification of unknown target genes of human transcription factors using chromatin immunoprecipitation.
Methods 2002, 26:3747. PubMed Abstract  Publisher Full Text

Hube F, Myal Y, Leygue E: The promoter competition assay (PCA): a new approach to identify motifs involved in the transcriptional activity of reporter genes.
Front Biosci 2006, 11:15771584. PubMed Abstract  Publisher Full Text

Kotelnikova EA, Makeev VJ, Gelfand MS: Evolution of transcription factor DNA binding sites.
Gene 2005, 347:255263. PubMed Abstract  Publisher Full Text

Dermitzakis ET, Clark AG: Evolution of transcription factor binding sites in Mammalian gene regulatory regions: conservation and turnover.
Mol Biol Evol 2002, 19(7):11141121. PubMed Abstract  Publisher Full Text

Dermitzakis ET, Bergman CM, Clark AG: Tracing the evolutionary history of Drosophila regulatory regions with models that identify transcription factor binding sites.
Mol Biol Evol 2003, 20(5):703714. PubMed Abstract  Publisher Full Text

Rajewsky N, Socci ND, Zapotocky M, Siggia ED: The evolution of DNA regulatory regions for proteogamma bacteria by interspecies comparisons.
Genome Res 2002, 12(2):298308. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Mustonen V, Lassig M: Evolutionary population genetics of promoters: predicting binding sites and functional phylogenies.
Proc Natl Acad Sci USA 2005, 102(44):1593641. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Gerasimova AV, Gelfand MS: Evolution of the NadR regulon in Enterobacteriaceae.
J Bioinform Comput Biol 2005, 3:10071019. PubMed Abstract  Publisher Full Text

Rodionov DA, Gelfand MS: Identification of a bacterial regulatory system for ribonucleotide reductases by phylogenetic profiling.
Trends Genet 2005, 21:385389. PubMed Abstract  Publisher Full Text

Rodionov DA, Vitreschak AG, Mironov AA, Gelfand MS: Comparative genomics of the regulation of methionine metabolism in Grampositive bacteria.
Nucleic Acids Res 2004, 32:33403353. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Gasch AP, Moses AM, Chiang DY, Fraser HB, Berardini M, Eisen MB: Conservation and evolution of cisregulatory systems in ascomycete fungi.
PLoS Biol 2004, 2(12):e398. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Costas J, Casares F, Vieira J: Turnover of binding sites for transcription factors involved in early Drosophila development.
Gene 2003, 310:21520. PubMed Abstract  Publisher Full Text

Mazon G, Campoy S, Erill I, Barbe J: Identification of the Acidobacterium capsulatum LexA box reveals a lateral acquisition of the Alphaproteobacteria lexA gene.
Microbiology 2006, 152(Pt 4):110918. PubMed Abstract  Publisher Full Text

Erill I, Jara M, Salvador N, Escribano M, Campoy S, Barbe J: Differences in LexA regulon structure among Proteobacteria through in vivo assisted comparative genomics.
Nucleic Acids Res 2004, 32(22):661726. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Rodionov DA, Dubchak IL, Arkin AP, Alm EJ, Gelfand MS: Dissimilatory metabolism of nitrogen oxides in bacteria: comparative reconstruction of transcriptional networks.
PLoS Computational Biology 2005, 1:e55. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Tsong AE, Tuch BB, Li H, Johnson AD: Evolution of alternative transcriptional circuits with identical logic.
Nature 2006, 443(7110):41520. PubMed Abstract  Publisher Full Text

Tanay A, Regev A, Shamir R: Conservation and evolvability in regulatory networks: the evolution of ribosomal regulation in yeast.
Proc Natl Acad Sci USA 2005, 102(20):72038. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Ihmels J, Bergmann S, GeramiNejad M, Yanai I, McClellan M, Berman J, Barkai N: Rewiring of the yeast transcriptional network through the evolution of motif usage.
Science 2005, 309(5736):93840. PubMed Abstract  Publisher Full Text

Rodionov DA, Gelfand MS, Todd JD, Curson ARJ, Johnston AWB: Computational reconstruction of iron and manganeseresponsive transcriptional networks in alphaproteobacteria.
PLoS Comput Biol 2006, 2(12):3163. Publisher Full Text

Permina EA, Gelfand MS: Heat Shock (sigma32 and HrcA/CIRCE) regulons in beta, gamma and epsilonproteobacteria.
J Mol Microbiol Biotechnol 2004, 6(34):174181. Publisher Full Text

Panina EM, Vitreschak AG, Mironov AA, Gelfand MS: Regulation of biosynthesis and transport of aromatic amino acid in lowGC Grampositive bacteria.
FEMS Microbiol Lett 2003, 222:211220. PubMed Abstract  Publisher Full Text

Gelfand MS: Evolution of transcriptional regulation networks in microbial genomes.
Curr Opin Struct Biol 2006, 16:420429. PubMed Abstract  Publisher Full Text

Nei M: Selectionism and neutralism in molecular evolution.
Mol Biol Evol 2005, 22(12):23182342. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Rockman MV, Hahn MW, Soranzo N, Goldstein DB, Wray GA: Positive selection on a humanspecific transcription factor binding site regulating IL4 expression.
Curr Biol 2003, 13(23):211823. PubMed Abstract  Publisher Full Text

Moses AM, Chiang DY, Kellis M, Lander ES, Eisen MB: Position specific variation in the rate of evolution in transcription factor binding sites.
BMC Evol Biol 2003, 3:19. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Kechris KJ, van Zwet E, Bickel PJ, Eisen MB: Detecting DNA regulatory motifs by incorporating positional trends in information content.
Genome Biol 2004, 5(7):R50. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Mazon G, Lucena JM, Campoy S, Fernandez de Henestrosa AR, Candau P, Barbe J: LexAbinding sequences in Grampositive and cyanobacteria are closely related.
Mol Genet Genomics 2004, 271:409. PubMed Abstract  Publisher Full Text

Mackiewicz P, ZakrzewskaCzerwinska J, Zawilak A, Dudek M, Cebrat S: Where does bacterial replication start? Rules for predicting the oriC region.
Nucleic Acids Res 2004, 32(13):37813791. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Rodionov DA, Mironov AA, Gelfand MS: Conservation of the biotin regulon and the BirA regulatory signal in Eubacteria and Archaea.
Genome Research 2002, 12:15071516. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Campoy S, Mazon G, Fernandez de Henestrosa AR, Llagostera M, Monteiro PB, Barbe J: A new regulatory DNA motif of the gamma subclass Proteobacteria: identification of the LexA protein binding site of the plant pathogen Xylella fastidiosa.
Microbiology 2002, 148(Pt 11):358397. PubMed Abstract  Publisher Full Text

Zverlov VV, Schwarz WH: Organization of the chromosomal region containing the genes lexA and topA in Thermotoga neapolitana. Primary structure of LexA reveals phylogenetic relevance.
Syst Appl Microbiol 1999, 22(2):1748. PubMed Abstract

Mazon G, Erill I, Campoy S, Cortes P, Forano E, Barbe J: Reconstruction of the evolutionary history of the LexAbinding sequence.
Microbiology 2004, 150(Pt 11):378395. PubMed Abstract  Publisher Full Text

Danilova LV, Gel'fand MS, Liubetskiĭ VA, Laĭkova ON: Computer analysis of regulating metabolism of glycerol3phosphate in proteobacteria genome.
Mol Biol (Mosk) 2003, 37(5):8439. PubMed Abstract

Jukes T, Cantor C: Evolution of protein molecules. In Mammalian Protein Metabolism. Academic Press; 1969:21132.

Kimura M: A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotides sequences.
J Mol Evol 1980, 16:111120. PubMed Abstract  Publisher Full Text

Galtier N, Gouy M: Inferring phylogenies from DNA sequences of unequal base compositions.
Proc Natl Acad Sci USA 1995, 92(24):1131721. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Galtier N, Gouy M: Inferring pattern and process: maximumlikelihood implementation of a nonhomogeneous model of DNA sequence evolution for phylogenetic analysis.
Mol Biol Evol 1998, 15(7):8719. PubMed Abstract  Publisher Full Text

Goncalves I, Robinson M, Perriere G, Mouchiroud D: JaDis: computing distances between nucleic acid sequences.
Bioinformatics 1999, 15(5):4245. PubMed Abstract  Publisher Full Text

Tamura K, Kumar S: Evolutionary distance estimation under heterogeneous substitution pattern among lineages.
Mol Biol Evol 2002, 19(10):172736. PubMed Abstract  Publisher Full Text

Gerland U, Hwa T: On the selection and evolution of regulatory DNA motifs.
J Mol Evol 2002, 55(4):386400. PubMed Abstract  Publisher Full Text

Gerland U, Moroz J, Hwa T: Physical constraints and functional characteristics of transcription factorDNA interaction.
Proc Natl Acad Sci USA 2002, 99(19):1201520. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Halpern A, Bruno W: Evolutionary distances for proteincoding sequences: modeling sitespecific residue frequencies.
Mol Biol Evol 1998, 15(7):9107. PubMed Abstract  Publisher Full Text

Kimura M: On the probability of fixation of mutant genes in a population.
Genetics 1962, 47:7139. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Kinney J, Tkacik G, Callan CJ: Precise physical models of proteinDNA interaction from highthroughput data.
Proc Natl Acad Sci U S A 2007, 104(2):5016. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Pollock DD, Bruno WJ: Assessing an unknown evolutionary process: effect of increasing sitespecific knowledge through taxon addition.
Mol Biol Evol 2000, 17:18541858. PubMed Abstract  Publisher Full Text

Perna NT, Kocher TD: Unequal base frequencies and estimation of substitution rates.

Collins TM, Wimberger PH, Naylor GJP: Compositional bias, characterstate bias, and characterstate reconstruction using parsimony.
Syst Biol 1994, 47:482496. Publisher Full Text

Keightley P, Kryukov G, Sunyaev S, Halligan D, Gaffney D: Evolutionary constraints in conserved nongenic sequences of mammals.
Genome Res 2005, 15(10):13738. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Lercher M, Smith N, EyreWalker A, Hurst L: The evolution of isochores: evidence from SNP frequency distributions.
Genetics 2002, 162(4):180510. PubMed Abstract  Publisher Full Text  PubMed Central Full Text