 Review
 Open Access
 Published:
From biophysics to evolutionary genetics: statistical aspects of gene regulation
BMC Bioinformatics volume 8, Article number: S7 (2007)
Abstract
This is an introductory review on how genes interact to produce biological functions. Transcriptional interactions involve the binding of proteins to regulatory DNA. Specific binding sites can be identified by genomic analysis, and these undergo a stochastic evolution process governed by selection, mutations, and genetic drift. We focus on the links between the biophysical function and the evolution of regulatory elements. In particular, we infer fitness landscapes of binding sites from genomic data, leading to a quantitative evolutionary picture of regulation.
Introduction
Genomic functions often cannot be understood at the level of single genes but require the study of gene networks. This systems biology credo is nearly commonplace by now. Evidence comes from the comparative analysis of entire genomes: Current estimates put, for example, the number of human genes at around 22000, hardly more than the 14000 of the fruit fly, and not even an order of magnitude higher than the 6000 of baker's yeast. The complexity and diversity of higher animals therefore cannot be explained in terms of their gene numbers. If, however, a biological function requires the concerted action of several genes, and conversely, a gene takes part in several functional contexts, an organism may be defined less by its individual genes but by their interactions. The emerging picture of the genome as a strongly interacting system with many degrees of freedom brings new challenges for experiment and theory, many of which are of a statistical nature. And indeed, this picture continues to make the subject attractive to a growing number of statistical physicists.
Genes encode proteins, and proteins perform functions in the cell. Hence, a gene takes part in a biological function only if it is expressed, i.e., if the protein produced from it is present in the cell. Genes interact by regulation: the protein of one gene can influence the production of protein from another gene. Gene regulation can take place during transcription, the process by which the cell reads the information contained in a gene and copies it to messenger RNA (which is subsequently used to make a functional protein). This is the most fundamental level of interactions between genes: the transcription of one gene may be enhanced or reduced by the expression of other genes. Transcriptional regulation is thus a good starting point for theory. We should keep in mind, however, that it is not the only mode of gene interactions. Especially in eukaryotes, additional regulation mechanims involving histones, chromatin, microRNAs etc. become relevant, which are just entering the stage of model building. An excellent introduction to the biology of regulation can be found in [1].
This article is a primer on theoretical aspects of gene interactions, and we limit ourselves to transcriptional regulation. Clearly, the subject has rather diverse aspects:
(1) Transcription is a biophysical process, which involves the interaction of DNA and proteins. Its regulation takes place through the binding of proteins to DNA at specific loci in the vicinity of the gene to be regulated. Already at this level, this process is rather complex and not yet fully understood. What enables the protein to find one or a few specific functional sites in a genome of up to billions of base pairs, bind there with sufficient strength to influence transcription, and leave again once its task is performed?
(2) Given that the protein can find its functional sites, can we as well? If that is possible, we can predict the specific gene interactions building regulatory networks from sequence data. The analysis of regulatory DNA is a major topic of research in bioinformatics, with the aim of identifying statistical characteristics of functional loci and of building search algorithms.
(3) Regulation is also becoming an important part of evolutionary biology [2, 3]. If regulatory networks are to explain the differentiation of higher animals, there must be efficient modes of evolution for the interactions between genes. At the level of regulary DNA, these modes remain largely to be explored. It is clear, however, that the underlying evolutionary dynamics is the basis of a quantitative understanding of regulatory networks.
All three aspects of regulation contribute to a unified theoretical picture. Key concepts such as the biophysical binding energy, the bioinformatic scoring function, and the evolutionary fitness turn out to be rather deeply related. We will focus on these crosslinks between different fields, which are likely to become important for future research. A challenge for an introductory presentation is the diversity of relevant background material, only a rather ecclectic account of which can be presented here. Yet, I hope it transpires even from this short introduction that present quantitative genomics is an area of science shaped by a remarkable confluence of ideas from different disciplines.
Biophysics of transcriptional regulation
The fundamental step in the regulatory interaction between two genes is a binding process: the protein produced by the first gene acts as a transcription factor for the second gene, i.e., it binds to a functional site on the DNA close to the second gene and thereby enhances or suppresses its transcription. Binding sites are short, typically segments of 10 to 15 base pairs in prokaryotes and even shorter segments in eukaryotes. They are primarily located in the cisregulatory region of a gene, which lies just upstream of its proteincoding sequence and extends over hundreds of base pairs in prokaryotes and over thousands of base pairs in eukaryotes. The scenario of transcriptional regulation is sketched in Fig. 1. A transcription factor bound to a functional binding site regulates the downstream gene by recruiting or repelling RNA polymerase. This proteinprotein interaction catalyzes or suppresses the process of transcription of the gene. All these binding processes should not be understood as on or off; they happen with certain probabilities, which are determined by the binding energies and the numbers of the molecules involved.
FactorDNA binding energies
The interaction of a transcription factor protein with DNA is twofold: There is a positionunspecific attraction with energy E_{ u }and a specific interaction, whose energy depends on the particular locus where the factor binds. The unspecific part is the electrostatic interaction between the positively charged protein and the negatively charged DNA backbone, while the specific part involves hydrogen bonds between the binding domain of the protein and the nucleotides of the binding locus. A locus is specified by its starting position r and its length ℓ (with relevant values ℓ of order 10). The specific binding energy E(r) depends on ℓ consecutive nucleotides a = (a_{1},..., a_{ℓ}) counted downstream from the starting position, the sequence state or genotype of that locus. Switching between unspecific and specific binding takes place via a conformation change of the factor protein. As a result of these interactions, the factor protein can be in three thermodynamic states as shown in fig. 2: unbound (i.e., freely diffusing), unspecifically bound (i.e., diffusing along the DNA backbone), and specifically bound.
The biophysics of factorDNA binding has been established in a series of seminal papers [4–7]. More recently, the characteristics of specific binding have been measured for some bacterial transcription factors [8–12]. These can be summarized as follows:

(a)
The single nucleotides of a binding locus a = (a_{1},..., a_{ℓ}) give approximately independent contributions to the binding energy,
$E(a)={\displaystyle \sum _{i=1}^{\ell}{\epsilon}_{i}({a}_{i}).}$

(b)
At each position i, there is typically one preferred nucleotide ${a}_{i}^{\ast}$ with ε_{ i }(${a}_{i}^{\ast}$) = min_{ a }ε_{ i }(a). Hence, there is a unique "ground state" sequence a* = (${a}_{1}^{\ast}$,...${a}_{\ell}^{\ast}$) with minimal binding energy E* = E(a*), i.e., with strongest binding.

(c)
Mismatches with respect to the minimumenergy sequence involve energy costs ε_{ i }(a)  ε_{ i }(${a}_{i}^{\ast}$) ≈ 1  3 k_{ B }T per nucleotide.

(d)
There is an energy difference E_{ u } E* ~15 k_{ B }T between unspecific and strongest specific binding. Experimental data for the binding energies ε_{ i }(a) are known only for a few transcription factors.
Approximate values for these energies can also be inferred from nucleotide frequencies in functional binding sites [10]. A promising recent approach is to infer binding energies from largethroughput expression data [13]. For orderofmagnitude estimates, one often uses the socalled twostate approximation [7], which is homogeneous in the nucleotide positions and distinguishes only between match and mismatch:
with ε ≈ 2k_{ B }T . In this approximation, the binding energy of a sequence a is simply related to the Hamming distance d(a, a*), i.e., the number of nucleotide mismatches between a and a*,
E(a) = E* + ε·d(a, a*).
Energy distribution in the genome
Fig. 3(a) shows the sequence of energy values E(r) found in a segment of the E. coli genome for a specific transcription factor, the cAMP response protein (CRP) This "energy landscape" looks quite random, i.e., consecutive energy values are approximately uncorrelated. The distribution W_{dat}(E) of energies over the entire noncoding part of the E. coli genome is shown in fig. 3(b). We can compare this with the distribution W_{0}(E) obtained from a random sequence with the same nucleotide frequencies (i.e., from a scrambled genome). According to eq. (1), the binding energy E is then a sum of independent random variables ε_{ i }, and its distribution becomes approximately Gaussian by the law of large numbers. Fig. 3(b) shows that the actual distribution W_{dat}(E) is indeed of the same form as W_{0}(E) for most energies. However, a closer look at the lowenergy tail of the distribution shows that there are significantly more strong binding sites than expected from a random sequence [14–16]. So at least some of them are there not by chance but for a reason.
Search kinetics
All three thermodynamic modes of a factor molecule – free diffusion, unspecific binding, and specific binding – are important for the search kinetics towards a functional site [4–6]. The unspecific attraction causes the transcription factor to be bound to DNA with a finite probability, i.e., a given molecule spends about equal amounts of time on and off the DNA backbone. Hence, the search process is a mixture of effectively onedimensional diffusion along the DNA backbone and threedimensional diffusion in the surrounding medium. This proves more efficient than purely one or threedimensional diffusion. In the 1D mode, the factor diffuses in a flat energy landscape if it is in the conformation of unspecific binding, or in the landscape E(r) if it is in the conformation of specific binding. In this way, it can sample the lowenergy part of the landscape E(r) while avoiding its barriers. The main obstacles on its way to a functional site are spurious binding sites, which have a low energy E(r) by chance and act as traps. We lack a completely satisfactory picture of the search kinetics, which is an area of current research [14, 17]. However, this process proves to be remarkably fast. Typical search times are less than a minute, i.e., substantially shorter than typical functional intervals in a cell cycle of at least minutes. Therefore, the regulatory effect of a site is related to its probability of binding a factor molecule at equilibrium, which can be evaluated by standard thermodynamics.
Thermodynamics of factor binding
We start with the idealized but instructive problem of a single factor protein interacting with a genome of length L ≫ 1, which contains a single functional site, while the rest of the sequence is random. Since the protein is bound to the DNA with a probability of about 1/2, we neglect the unbound state for the subsequent probability estimates and study only the bound protein, which is at equilibrium between specific and unspecific binding. At each position r, the likelihood of these two states is given by the Boltzmann factors exp [E(r)/k_{ B }T] and exp [E_{ u }/k_{ B }T], respectively. Hence, the partition function for a single protein has the form
The functional site, which is assumed to be positioned at r = r_{ f }, must have a low specific binding energy E = E(r_{ f }). We now single out this position and write
where Z_{0} is the partition function of a completely random sequence. The probability of the factor being bound specifically at the functional site is then
where F_{0} = k_{ B }T log Z_{0} is the free energy for a random genome. Thus, the binding probability depends on the binding energy in a sigmoid way, with a threshold energy E = F_{0} between strong and weak binding.
This strongly nonlinear dependence is known to physicists as a Fermi function.
It is easy to generalize the thermodynamic formalism to more than one factor molecule. Ignoring the overlap between close sites, each position r can be empty or be occupied either by an unspecifically or by a specifically bound factor. Using a chemical potential σ, the manyfactor partition function can hence be written as
where
is a sum over the three thermodynamic states at position r: no factor bound, one factor bound specifically or unspecifically. The chemical potential σ is determined by the number of factor molecules, n, via the relation n = (d/dσ) log Z(σ). For actual transcription factor numbers, which are of order 1  10^{4}, this relation is well approximated by [14]
The functional site is now occupied by a specifically bound factor with probability
The binding probability – and hence the effects of the functional site on the regulated gene – are thus determined by the binding energy, the number of factor molecules, and on the genomic background (via the free energy F_{0}). The dependence p(E) is a Fermi function with threshold energy E = F_{0} + k_{ B }T log n, which is shifted with respect to the singlemolecule case. Clearly, p is also a Fermi function of log n at fixed binding energy, with a threshold at log n = (E  F_{0})/k_{ B }T . If there is more than one functional site in the genome, the calculation remains unaffected as long as their number is much smaller than n.
Sensitivity and genomic design of regulation
The regulatory machinery can be very efficient: in bacteria, it has been shown that single factor molecules can have regulatory effects. We can use eq. (6) to enquire how the cell can reach this high level of sensitivity, following mostly ref. [14]. We assume a minimal genome which has a single functional site of maximum binding strength E* and is otherwise random. If a single factor molecule is to affect regulation, its binding to the functional site must not be overwhelmed by the remainder of the genome. This leads to a criterion on the signaltonoise ratio of regulatory interactions,
F_{0} ≳ E*, (11)
which in turn imposes a number of constraints on the design of regulatory DNA:

(a)
In a random genome, there must be at most a number of order one minimumenergy binding sites. Estimating the probability to find such a site at a given position as (1/4)^{ℓ}, we obtain the condition
L(1/4)^{ℓ} ≲ 1. (12)
This gives a lower bound on the site length, ℓ ≳ log L/log 4. For a bacterial genome (L ~10^{6}), we obtain ℓ ≳ 10, which gives the right length of functional binding sites. However, this bound is not fulfilled in eukaryotes. Indeed, eukaryotic genomes use a different design with groups of adjacent binding sites.

(b)
For each minimumenergy site, there are ℓ suboptimal sites of Hamming distance 1 from the minimumenergy sequence. These must not suppress the binding to the minimumenergy site, i.e.,
exp(E*/k_{ B }T) ≳ ℓ exp [(E* + ε)/k_{ B }T] (13)
in the twostate approximation. This gives a lower bound on the binding energy per nucleotide, ε/k_{ B }T ≳ log ℓ ≈ 2  3.

(c)
Finally, the unspecific binding in the entire genome must not suppress the specific binding to a minimumenergy site, i.e.,
exp(E*/k_{ B }T) ≳ L exp(E_{ u }/k_{ B }T). (14)
This produces a lower bound on the energy gap between unspecific and optimal specific binding, (E_{ u } E*)/k_{ B }T ≳ log L ≈ 15.
Quite remarkably, these bounds are fulfilled as approximate equalities in bacteria. Hence, the machinery of transcriptional regulation operates just at the treshold of singlemolecule sensitivity, i.e, F_{0} ≈ E*.
Programmability and evolvability of regulatory networks
Of course, not every regulatory interaction is equally sensitive. To switch genes on or off, the cell uses the dependencies of the binding probability both on factor numbers and on binding energies. During the cell cycle, the level of n can vary over several orders of magnitude, say, between a few and tens of thousands of molecules. At a given value of n, the effects on the regulated genes differ since their functional sites have different values of E. The binding energies can change on evolutionary time scales by mutations of the site sequence, which leads to regulatory differences between individuals and, ultimately, between species. Both parameters are thus necessary to encode pathways in regulatory networks. This is most flexible if minimumenergy sites are indeed sensitive to a single factor molecule as discussed above. Differential programmability as a network design principle [14] thus favors complicated molecular structures with longer binding sites and larger binding energies. However, this competes with the evolvability of the system by a stochastic evolution process [18]. We have seen that the singlemolecule sensitivity is just marginally reached in bacteria. This indicates that the actual machinery may result from a compromise between programmability and evolvability: binding sites are just complicated enough to work. It also indicates that genomic structures can only be understood from their evolution. This aspect will be developed further below, after we have introduced sequence analysis aspects in the next section.
Bioinformatics of regulatory DNA
Predicting regulatory interactions between genes is clearly a key problem in bioinformatics, which is as important as the analysis of individual genes and proteins. It is not surprising that this problem is very difficult since, as we have discussed in the previous section, targeting regulatory input in a large genome is a tremendous signaltonoise problem even for the cell itself. Its solution via the analysis of regulatory DNA requires finding statistical criteria to distinguish between functional binding sites and background sequence. A general introduction to the relevant sequence statistics can be found in ref. [19].
Markov model for background sequence
We begin by specifying a stochastic model for the nonfunctional segments of intergenic DNA. These are assumed to be Markov sequences with uniform singlenucleotide frequencies p_{0}(a) (a = A, C, G, T). Hence, the probability of finding a given sequence has the factorized form
This assumption should not be taken too literally. The term "nonfunctional" refers to binding of a particular transcription factor. Intergenic DNA contains plenty of nonrandom elements with other functions (e.g., binding sites for other factors) or without known function (such as repeat elements). The salient point is, however, that most of intergenic DNA is well approximated by a Markov sequence with respect to binding of a given transcription factor. To make this more precise, we project the distribution P_{0}(a) for segments of length ℓ onto the binding energy E as independent variable. Denoting the projected distribution for simplicity with the same letter P_{0}, we have
This distribution is close to the actual genomic distribution W_{dat}(E) for most values of E, as we have seen in fig. 3. It is possible to improve the background model by introducing small frequency couplings between neigboring letters [15, 16].
Probabilistic model for functional sites
The sequences a = (a_{1},..., a_{ℓ}) at functional sites of a given transcription factor are assumed to be drawn from a different distribution Q(a). We write this distribution in the form
Q(a) = P_{0}(a) exp [S(a)]. (17)
The quantity S(a), which is called the relative log likelihood score of the distributions P_{0} and Q, will turn out to have an important evolutionary meaning as well.
The singlenucleotide distribution q_{ i }(a) at a given position i within functional loci is obtained by summing the full distribution Q over all other positions
The set of these marginal distributions, q_{ i }(a) (i = 1,..., ℓ; a = A, C, G, T) is called the position weight matrix for binding sites of a given factor [20]. If the score function is additive in the nucleotide positions, $S(a)={\displaystyle {\sum}_{i=1}^{\ell}{s}_{i}({a}_{i})}$, the Q distribution has a factorized form, $Q(a)={\displaystyle {\prod}_{i=1}^{\ell}{q}_{i}({a}_{i})}$with
q_{ i }(a) = p_{0}(a) exp [s_{ i }(a)]. (19)
This additivity assumption is made in most of the existing literature since the position weight matrix (18) can be inferred from a sample of known functional site sequences, which in turn determines directly the single nucleotide scores (19). This scoring is the basis for a number of site prediction methods in single species and by crossspecies analysis; see, e.g., refs. [20–24].
Here we treat functional sites as coherent statistical units and do not make the assumption of additivity of the score function [16]. As will be discussed in the next section, functionality imposes correlations between the nucleotide frequencies within a functional site, preventing factorization of the Q distribution. Of course, it is not possible to reconstruct the full distribution Q(a), which lives on a 4^{ℓ}dimensional sequence space, from a limited sample of experimentally known functional sites. However, we can again project this distribution onto the binding energy as independent variable, Q(E) = ∑_{ a }Q(a)δ(E  E(a)). Since all regulatory effects of a functional site depend on its sequence a only via the binding energy, we can also write the score as a function of the energy, S(a) = S(E(a)) (this will become obvious in the next section). Hence, the relationship (17) has the same form for the projected distributions,
Q(E) = P_{0}(E) exp [S(E)]. (20)
Bayesian model for genomic loci
Assuming that functional loci are distributed randomly with a small probability λ, we now combine the models for background sequence and for functional sites into a model for the full distribution of sequences a in intergenic DNA,
W(a) = (1  λ) P_{0}(a) + λ Q(a). (21)
(At the moment, we are ignoring the possible overlap between functional sites). In the language of statistics, this is a probabilistic model with hidden variables. The output of this model consists of pairs (m, a): First, the model variable m is randomly drawn, labelling a locus as nonfunctional (m = 0) with probability 1  λ or as functional (m = 1) with probability λ. Then the sequence is drawn from the corresponding distribution P_{0}(a) or Q(a). However, only the sequence counts a are available data. The "hidden" variable m can be inferred from the data in a probabilistic way using Bayes' formula, which expresses the joint probability distribution of data and model in terms of its conditional and its marginal distributions
prob(a, m) = prob(am) prob(m) = prob(ma) prob(a) (23)
with prob(a) = ∑_{ m }prob(am)prob(m). We can solve for the conditional probability of the model for given data a,
For the probability of functionality, ρ_{ f }(a) = prob(m = 1a), this formula reads
The dependence on S has again the form of a Fermi function. Its threshold value S = log [(1  λ)/λ] separates sequences that are more likely to be functional or more likely to be background.
The full Bayesian model (21) can again be projected onto the energy variable,
W(E) = (1  λ)P_{0}(E) + λ Q(E).
In this form, it can be tested against genomic data [16]. To plot the distributions P_{0}, Q, and W as functions of E, we use eq. (1) with an energy matrix ε_{ i }(a) = ε_{0} log [q_{ i }(a)/p_{0}(a)] estimated from the position weight matrix up to an overall constant ε_{0} [10]. For our example of the CRP transcription factor, the distribution Q(E) can be estimated from the about 50 known binding sites in the E. coli genome. Using this Q distribution and a probability of functionality λ ≈ 6 × 10^{4}, the full distribution W(E) produces an excellent fit of the count histogram W_{dat}(E) over the entire range of energies; see fig. 4(a). The log likelihood score function S(E) = log [Q(E)/P_{0}(E)] is shown in fig. 4(b), shifted such that the curve has its zero at a point E_{ s }≈ 13 beyond which binding becomes negligible.
The resulting probability of functionality ρf(E) as given by eq. (24) is also shown in fig. 4(b). This indicates the dilemma for the prediction of individual binding sites based on sequence data from a single species. Many functional sites have energies in the "twilight" region between the ensembles λQ and (1  λ)P_{0}, where ρf takes values around 1/2. Hence, depending on the energy cutoff chosen, any prediction is torn between many false negatives or many false positives.
Dynamic programming and sequence analysis
It is straightforward to generalize the Bayesian approach to longer segments of intergenic DNA, which are covered by an unknown number s of nonoverlapping functional sites as shown in fig. 5[22]. The hidden variables are now the sequence of left initial positions r_{ f }≡ (r_{1},..., r_{ s }) of the functional sites (with the nooverlap constraint r_{ν+1}≥ r_{ ν }+ ℓ for ν = 1,..., s  1). The full sequence distribution in a segment of length L has the form
where Z is a normalization factor, $\tilde{\lambda}$ = λ + O(λ^{2}) is a weight factor for each functional locus (the negligible correction terms originate from the nooverlap constraint), and W_{ L }(a_{1},..., a_{ L }r_{ f }) is the sequence distribution for given positions of functional loci,
with r_{n+1}= L + 1. The sum over sequences r_{ f }of arbitrary length s seems formidable at first, but W_{ L }is easy to compute from the recursion
W_{ r }(a_{1},..., a_{ r }) = (1  $\widehat{\lambda}$)p_{0}(a_{ r })W_{r1}(a_{1},..., a_{r1}) + $\tilde{\lambda}$Q(a_{rℓ+1},..., a_{ r }) W_{rℓ}(a_{1},..., a_{rℓ}) (28)
with the initial condition W_{0} = 1 and $\widehat{\lambda}$ = $\tilde{\lambda}$ + O($\tilde{\lambda}$^{2}). This type of recursion relation is usually called a dynamic programming algorithm in computer science. In physics, it is known as a transfer matrix, and the sum (27) is recognized as the corresponding discrete path integral in imaginary time r, if we interpret r_{ f }as encoding a path m(r) that takes the value m = 1 at the nucleotide positions r_{ ν },..., r_{ ν }+ ℓ  1 (ν = 1,..., s) within functional loci and m = 0 otherwise (see fig. 5). Both concepts prove very useful also in more general problems of sequence alignment.
In analogy to (24), the probability of a set r_{ f }of functional loci for given sequence data is
The most likely set ${r}_{f}^{\ast}$ can be obtained by the following "backward" algorithm: Given the sequence (W_{1},..., W_{ L }) obtained from the "forward" recursion (28), we can decide for every point r whether it is more likely to be a background position or the endpoint of a functional locus, ignoring all sequence information from positions > r. This depends on whether the leading contribution to W_{ r }comes from the first or second term on the r.h.s. of (28) and defines the local optimum model m*(r). The global optimum set of functional loci respecting the nooverlap constraint is then ${r}_{f}^{\ast}$ = {rb(r) = 1}, where b(r) is given by the recursion b(r) = ℓ if b(r + 1) ≤ 1 &m*(r) = 1 and b(r) = max(b(r + 1)  1, 0) otherwise, with the initial condition b(L + 1) = 0.
The Bayesian model can easily be extended to sequences containing several types of binding sites, which bind different transcription factors and are distinguished by their Q distributions. Dynamic programming algorithms can thus predict the likely coverage of a sequence with binding sites of known type [22]. This is the first step in extending the statistical analysis from single binding sites to entire regions of regulatory DNA. Indeed, models of this kind have been applied successfully to predict regulatory elements in eukaryotes, which typically consist of functional groups of adjacent binding sites. In the algorithms currently used, however, the scoring in (27) is strictly additive for groups of nonoverlapping binding sites: it does not take into account dependencies between the sites within one functional group or overlapping sites within one sequence.
Evolution of regulatory DNA
In the statistical picture developed so far, background sequences and functional sites are reduced to ensembles P_{0} and Q. This picture is incomplete in two ways. On one hand, it is quite disconnected from the biophysical aspects discussed before: the specific function of binding sites hardly enters the standard formalism of position weight matrices. On the other hand, there is not yet any notion of time and dynamics. Sequences change by various mutation processes, and the observed sequence ensembles derive from this evolutionary dynamics. The evolution of functional loci is fundamentally different from that of background sequence: it is subject to natural selection, that is, the fitness of an organism depends on its genotype a at a functional locus via the effects on the regulated gene. At this point, the biophysics of binding enters the evolution of functional sequences [25–27]. Moreover, it becomes clear that the statistical framework has to be extended from individual sequences to distributions of genotypes in a population. In this section, we develop an evolutionary picture of regulatory DNA, from which we obtain expressions for the sequence ensembles P_{0}, Q, and the score function S. The next four paragraphs are a selfcontained introduction to the underlying concepts of population genetics.
Deterministic population dynamics and fitness
We start by describing the evolution of a large population, which contains individuals of different genotypes a. Each genotype is assumed to produce a specific phenotype, which may influence the reproductive success of the individuals carrying it. With respect to factor binding, the phenotype can be associated with the binding energy E(a), since presumably all organismic effects of a locus depend on its genotype only via the binding energy. However, the discussion in the following paragraphs is more general.
We first assume that the subpopulations of a given genotype reproduce separately, i.e., there neither transitions between genotypes through mutations nor (in a sexually reproducing population) mixing through genomic recombination. Writing the dynamics of the subpopulations in the form of simple growth laws,
defines the (Malthusian) fitness F_{ a }(t) of each genotype. For notational simplicity, we now limit ourselves to the case of just two genotypes a and b, where (30) can be written as growth laws for the total population size N(t) = N_{ a }(t) + N_{ b }(t) and for the population fraction x(t) = N_{ b }(t)/N (t) of genotype b,
with $\overline{F}$(t) ≡ [1  x(t)]F_{ a }(t) + x(t)F_{ b }(t) and ΔF_{ ab }(t) ≡ F_{ b }(t)  F_{ a }(t). This decomposition is useful since the overall growth rate $\overline{F}$(t) is often strongly timedependent due to external conditions (e.g., seasonality), while fitness differences, which reflect intrinsic properties of the phenotypes, are more stable. Different genotypes coexisting in a population frequently produce the same or very similar phenotypes and thus have equal fitness (ΔF_{ ab }= 0).
Assuming ΔF_{ ab }to be constant over the time of observation, the solution of eq. (32) is the evolutionary trajectory
with the initial condition x(t_{0}) = x_{0}, shown in fig. 6(a). For ΔF_{ ab }≠ 0, the fixed points of this dynamics are the monomorphic population states x = 0, and x = 1, of which x = 1 is stable for ΔF_{ ab }> 1 and x = 0 for ΔF_{ ab }< 1. The approach to the stationary state takes place on a characteristic time scale τ_{ d }= 1/ΔF_{ ab }.
In the important case of neutral evolution (ΔF_{ ab }= 0), the evolutionary outcome remains indefinite. These results, which can readily be generalized to more than two phenotypes, are a simple version of Fisher's fundamental theorem of natural selection: any population with initially coexisting phenotypes of different fitness will evolve towards a state where only the fittest phenotype is present.
Fisher's theorem seems to prove the popularized Darwinian notion of the "survival of the fittest". However, it rests on very restrictive assumptions that are never fulfilled in a natural population. The deterministic growth law (32) neglects mutations and recombinations, as well as the reproductive fluctuations present in any population due to its finite number of individuals. These other evolutionary forces have to be incorporated in our theoretical picture before we can even define fitness as a measurable quantity and before the theory can address the important case of neutral evolution.
Stochastic dynamics and genetic drift
Stochastic fluctuations of the reproduction process in a large but finite population have been studied extensively in population genetics, see [28, 29]. They are called genetic drift, an unfortunate name which may falsely suggest a deterministic effect. To take these fluctuations into account, we replace eq. (30) by a stochastic growth law,
where χ_{ a }(t) are Gaussian random variables with $\overline{{\chi}_{a}(t)}=0$ and
This form of noise is simply due to the law of large numbers, and the continuum dynamics (34) emerges as an effective largeN description for a plethora of discrete evolution models, which are defined at the level of individuals and have finite generation times. In the application to real populations, N has to be interpreted as the socalled effective population size, which can be inferred from genome data and is in general smaller than the actual population size.
In the case of two genotypes, eq. (34) can again be projected onto the population fraction x,
where χ_{ x }(t) = (∂x/∂N_{ a })χ_{ a }(t) + (∂x/∂N_{ b })χ_{ b }(t) are Gaussian random variables with zero mean and
This dynamics produces stochastic evolutionary trajectories x(t) as shown in fig. 6(b). To capture their statistics, we convert the Langevin equation (36) into a FokkerPlanck equation for the probability distribution of the genotype composition [28, 30],
The mathematical subtlety of this equation lies in the xdependent diffusion "constant" x(1  x)/2N, which reflects the multiplicative nature of the reproduction process. As a consequence, the two monomorphic population states x = 0 and x = 1 are also fixed points also of the stochastic dynamics. Any evolutionary trajectory x(t) will eventually lead to one of these states with probability 1; this is called the fixation of the corresponding genotype in the population. In other words, the FokkerPlanck equation (38) describes diffusion in the interval (0, 1) with absorbing boundaries. There is a family of stationary states
(x) = (1  φ)δ(x) + φδ (1  x),
parametrized by the fixation probability φ of genotype b. The value of φ depends on the initial condition x_{0} and can be computed by solving the backward diffusion equation
For timeindependent ΔF_{ ab }, the stationary solution φ (x_{0}) = lim_{t→∞}$\mathcal{\text{P}}$(x = 1, tx_{0}, t_{0}) has the form [28, 30]
which for nearneutral evolution (NΔF_{ ab }≪ 1) reduces to
φ(x_{0}, 0, N) = x_{0} + NΔF_{ ab }x_{0}(1  x_{0}) + ....
The characteristic time τ_{ s }of the stochastic dynamics interpolates between the diffusive scale N and the deterministic scale: τ_{ s }≈ min(N, τ_{ d }). It determines the typical time of the evolution process up to fixation, shown shaded in fig. 6(b).
Hence, the stochastic population dynamics depends no longer only on the fitness difference of the genotypes as in the deterministic case, but also on the initial state of the population and the the population size. Yet, our evolutionary picture is still incomplete. Population states with coexisting genotypes enter the dynamics as initial conditions, but since mutations are neglected, the model does not explain how this coexistence is generated and maintained.
Mutation processes and evolutionary equilibria
At the level of an individual, mutations are rare stochastic genotype changes a → b, which take place with rates μ_{a→b}, often coupled to the reproduction process. (These rates are all of the same order of magnitude, in estimates we therefore omit the indices.) We include mutations into the population dynamics (34) by their systematic effect on the genotype subpopulations,
while their stochastic effect (whose variance is of order Nμ) is neglected since it is small against the reproductive sampling noise χ_{ a }(t). In the case of two different genotypes, this dynamics can again be projected onto the variable x,
which leads to the FokkerPlanck equation [31]
For timeindependent ΔF_{ ab }, this equation has a single stable stationary state,
with a normalization constant Z that can be expressed in terms of Bessel and Gamma functions [32].
Substitution dynamics
Here we are interested in the stochastic evolution (45) and its equilibrium state (46) for Nμ ≪ 1, which is the relevant dynamical regime for nuclear DNA in eukaryotes and in most prokaryotes (but not in viral systems). In this regime, the mutation term in (45) is small against the diffusion term except for values of x close to the boundaries 0 or 1. In this region, the continuum approximation of eq. (45) is no longer valid, and (46) has to be replaced by a stationary solution $\mathcal{\text{P}}$_{ d }(N_{ a }) of the underlying discrete evolution model, which gives the probability that the population contains N_{ a }individuals of genotype a (with N_{ a }= N  N_{b} = 0, 1,..., N). The discrete solution is easily shown to have the singularity ${\mathcal{\text{P}}}_{d}(0)\simeq {(N{\mu}_{a\to b})}^{1}{P}_{d}(1)$. This singularity is correctly captured if we use the approximation ${P}_{d}({N}_{a})\simeq {\displaystyle {\int}_{{N}_{a}/N}^{({N}_{a}+1)/N}dx\mathcal{\text{P}}(x)}$ for all N_{ a }(except at the other boundary, where there is a similar singularity ${\mathcal{\text{P}}}_{d}(N)\simeq {(N{\mu}_{b\to a})}^{1}{P}_{d}(N1))$[33].
From this solution, we read off the following characteristics of the evolutionary dynamics at equilibrium, which are illustrated by the trajectory of fig. 6(c)[32]:

(a)
For sufficiently small values of μ, the population remains monomorphic for most of the time. Using the shorthands Q(a) = $\mathcal{\text{P}}$_{ d }(N_{ a }= 0) and Q(b) = $\mathcal{\text{P}}$_{ d }(N_{ a }= N), we have
Q(a) + Q(b) = 1  O(μN log N).

(b)
The ratio of probabilities for the two monomorphic population states is given by the ratio of "forward" and "backward" mutation rate, the fitness difference, and the effective population size:
$$\frac{Q(b)}{Q(a)}=\frac{{\mu}_{a\to b}}{{\mu}_{b\to a}}\mathrm{exp}(2N\Delta {F}_{ab})+O({N}_{\mu}).$$(48) 
(c)
The monomorphic population states x = 0 and x = 1 are unstable due to mutations even at arbitrarily small values of μ, which cause occasional transitions of the entire population from genotype a to b, and vice versa. These socalled substitutions are marked by dashed lines in fig. 6(c). The substitution rate u_{a→b}can be evaluated as the product of creating a single mutant of genotype b in an initially monomorphic a population, $N{\mu}_{a\to b}$, and its probability of fixation, φ (x_{0} = 1/N, ΔF_{ ab }, N). The time between initial mutation and fixation (shown by grey shading in fig. 6(c)) is still of order τ_{ s }and thus much shorter than the time scale 1/μ, on which mutation effects become important. Hence, the fixation probability φ is given to leading order by (41), which has been derived for μ = 0. Together we have [28, 30]
$${u}_{a\to b}=N{\mu}_{a\to b}\frac{1\mathrm{exp}(2\Delta {F}_{ab})}{1\mathrm{exp}(2N\Delta {F}_{ab})}.$$(49)
Hence, the substitution rate u_{a→b}is enhanced over μ_{a→b}for ΔF_{ ab }> 0 and suppressed for ΔF_{ ab }< 0, as shown in Fig. 7. For weak selection (NΔF_{ ab } ≪ 1), eq. (49) becomes
u_{a→b}= μ_{a→b}(1 + NΔF_{ ab }+ ...).
This reproduces Kimura's famous original result: for neutral evolution, the substitution rate equals the mutation rate in an individual, independently of the population size. For this reason, the rates μ_{a→b}are referred to as neutral mutation rates. For strong selection (NΔF_{ ab } ≫ 1 ≫ ΔF_{ ab }), eq. (49) takes the asymptotic forms
The backward substitution rate u_{b→a}is given by a formula similar to (49) with ΔF_{ ba }= ΔF_{ ab }. Forward and backward substitution rate have the simple ratio
for N ≫ 1. Comparing with (48), we obtain the consistency condition
Hence, for sufficiently small mutation rates (μN log N ≪ 1), a simple picture emerges: The evolution of a population can be described as a sequence of transitions between monomorphic genotype states (substitutions). The substitution rate u is determined by the corresponding mutation rate in an individual, the fitness difference between the genotypes, and the effective population size.
Neutral dynamics in sequence space, sequence entropy
This evolutionary picture can be generalized to multiple genotypes, for example, the 4^{ℓ} dimensional sequence space of genomic loci a = (a_{1},..., a_{ℓ}). Transitions between different sequence states are point mutations a → b, which change exactly one nucleotide. (We neglect here insertion and deletion processes, which change the length of the sequence). We first discuss neutral evolution, where the substitution rate u_{a→b}equals the mutation rate in an individual, μ_{a→b}, for all elementary transitions a → b. Bona fide neutral mutation rates can be inferred from DNA sequence alignments of sufficiently close species, recent insights have also come from studying repeat elements.
We assume the neutral dynamics has an equilibrium distribution P_{0}(a) which obeys detailed balance, i.e., the relation
holds for each pair of sequence states linked by an elementary transition process a → b. This says that the probability current at equilibrium, μ_{a→b}P_{0}(a)  μ_{b→a}P_{0}(b), vanishes for each elementary transition. Clearly, any distribution P_{0}(a) satisfying the conditions (54) is stationary under the dynamics with rates μ_{a→b}, but not every such dynamics has a stationary distribution which satisfies (54) (the simplest counterexample involving three states and a circular probability current a → b → c at stationarity).
However, as will be verified below, detailed balance is a good approximation for the genomic substitution dynamics at least in prokaryotes. (There are known violations at CpG islands in eukaryotes [34]). In the simplest type of models, every nucleotide a mutates independently of all other positions with uniform rates μ_{a→b}(i.e., μ_{a→b}= μ_{a→b}for any two sequences a = (..., a ,...) and b = (..., b, ...) differing by exactly one nucleotide). This produces a factorized equilibrium distribution P_{0}(a) of the form (15).
We can project the equilibrium distribution onto a measurable quantity as independent variable. For binding site sequences, a convenient choice is the binding energy E, and the projected distribution P_{0}(E) has the form (16). Hence we can define the sequence entropy [35]
S_{0}(E) = log P_{0}(E), (55)
which counts the log density of sequence states a at energy E, weighed by the distribution P_{0}(a).
Dynamics under selection, the scorefitness relation
The dynamics of substitutions can be studied in the same way for evolution under selection, which is specified at the level of genotypes by an arbitrary fitness function F(a) [18, 36]. This generalizes the results of [37] for a model with selection acting independently at different nucleotide positions, i.e., $F(a)={\displaystyle {\sum}_{i=1}^{\ell}{f}_{i}({a}_{i})}$. For each elementary transition a → b, the substitution rate u_{a→b}is determined by the neutral rate μ_{a→b}, the fitness difference ΔF_{ ab }, and the effective population size N according to (49). Given the detailed balance (54) of neutral evolution and the relation (52) between forward and backward rates, it then follows immediately that the evolutionary dynamics under selection also obeys detailed balance, as given by (53) with an equilibrium distribution Q(a) of the form (48). Thus we have [18, 36]:
The equilibrium distribution Q(a) of fixed genotypes generated by a substitution dynamics (49) with fitness function F(a) is related to its neutral counterpart P_{0}(a) by
Q(a) = P_{0}(a) exp [2NF (a) + const.], (56)
with the constant given by normalization.
We can project eq. (56) onto the fitness as independent variable. Defining the distribution Q(F) ≡ ∑_{ a }Q(a)δ (F (a)  F), similarly P_{0}(F), and the sequence entropy S_{0}(F) = log P_{0}(F), the projected identity takes the form
Q(F) = exp [2NF + S_{0}(F) + const.] (57)
For binding site sequences, we have a similar projection on the binding energy, Q(E) = exp [2NF(E) + S_{0}(E) + const.], since all genotypes with the same "phenotype" E have the same fitness, i.e., the same score S. The projected identities express the equilibrium distribution under selection in terms of fitness and sequence entropy, reflecting the balance between stochasticity (genetic drift) and selection [18]. For strong selection, the exponent 2NF  S_{0} is dominated by the fitness term, and Q(F) takes appreciable values only at points of nearmaximal fitness, i.e., where F_{ max } F ≲ 1/2N. For moderate selection, there is a nontrivial balance between both terms, and for weak selection, the Q distribution can be approximated by its neutral counterpart P_{0} = exp(S_{0}). Clearly, the roles of fitness and sequence entropy are formally analogous to those of energy and entropy in statistical physics of thermodynamic systems, if 2N is identified with the inverse temperature 1/k_{ B }T . Some consequences of this analogy are discussed in ref. [38].
The dynamics of substitutions establishes a rather general evolutionary grounding of genome statistics, if we identify the equilibrium distributions P_{0}(a) and Q(a) with the genomic distributions discussed in the previous section, as already anticipated by our notation. Comparing eqs. (56) and (17) gives a relation between fitness and score [16, 18]:
The loglikelihood score S(a) = log [Q(a)/P_{0}(a)] equals the fitness function multiplied by twice the effective population size up to a constant,
S(a) = 2NF(a) + const.. (58)
This relation allows us to use sequence data of a given genome to infer quantitative patterns of its evolution. We now discuss specific consequences for the evolution of regulatory DNA; an application to protein evolution can be found in ref. [37].
Measuring selection for binding sites
We first give a precise definition of functionality for regulatory (and other) elements: A binding locus is functional if the genotype at that locus is under selection (for binding of the corresponding factor). Nonfunctional loci have evolutionarily neutral genotypes. This definition asks whether binding at a given locus makes a difference to the organism or not. It is weaker than that of a functional binding site, which is a functional locus with a sequence a that is likely to actually bind the factor. A functional locus can lose its binding sequence due to deleterious mutations, leading to suboptimal fitness of the organism. Conversely, a nonfunctional locus can have by chance a sequence which does bind the factor: this is a spurious binding site without consequences for the organism. To measure the selection on functional sites in silico, we apply the identity (58) to the genomic distributions P_{0}(a) and Q(a). (Assuming equilibrium for most loci seems to be justified for our example of CRP binding sites in E. coli since we find very similar distributions in the distant bacterial species Salmonella typhimurium, and the factor protein itself is highly conserved between these species.) After projection onto the energy, the fitness landscape 2NF(E) for CRP binding sites is thus given by fig. 4(b)[16]. The fitness is constant in the nobinding region (E ≳ E_{ s }≈ 13) since the evolution is always neutral in that region. This constant is set to 0 in our normalization, i.e., F(E) measures the fitness gain of functional sites due to factor binding. Loci with strong binding are also under strong selection, with effective fitness values 2NF of order 10. Genetic drift counteracts selection, producing also loci with weaker binding and reduced effective fitness. This fitness "landscape" is thus qualitatively of the form predicted from the underlying biophysics [18, 25]. Of course, it should be kept in mind that this landscape results from averaging over a family of binding sites, which may have a spectrum of individual selection coefficients and selected binding strengths.
Nucleotide frequency correlations
A further consequence of (57) is the generic occurence of nucleotide frequency correlations within functional loci [18]. If the fitness function F(a) is not additive in the nucleotide positions, nucleotide frequencies are correlated in selected genotypes even if they are independent under neutral evolution. This happens quite generically since selection acts on the entire genotype a as a functional unit and not on its single nucleotides. For binding sites, fitness effects follow from the expression level of the regulated gene, which depends on the sequence a via the binding probability of the corresponding transcription factor. While the binding energy is often approximately additive in the nucleotide positions as given by (1), the binding probability (10) is a strongly nonlinear function of the energy. This introduces correlations between nucleotide frequencies at any two positions within functional loci, preventing factorization of the distribution Q(a).
Stationary evolution of binding sites
Functional loci with a substantial level of selection (as found for the CRP binding sites in E. coli) evolve in a way quite different from background sequence. This is quantified in fig. 8(a), which shows pairs of binding energies (E_{1}, E_{2}) for experimentally verified CRP binding sites in E. coli and the corresponding sites regulating orthologous genes in S. typhimurium [16, 27]. The evolutionary distance t between the two species and characteristics of the neutral mutation process can be inferred from alignments of background sequence. The "phenotypic" evolution of CRP binding is quantified by the energy transition probabilities G_{0}(E_{2}E_{1}) under neutral evolution and G_{ f }(E_{2}E_{1}) under stationary selection [16]. These are readily obtained by simulating the substitution dynamics over a time interval t for given initial value E_{1}, both with neutral rates μ_{a→b}and with rates u_{a→b}given by (49) and the fitness function 2NF(E) measured in E. coli. The resulting conditional expectation values ⟨G_{0}(E_{2}E_{1})⟩ and ⟨G_{ f }(E_{2}E_{1})⟩ for the binding energy in S. typhimurium are also shown in fig. 8(a). The data conform to the selection model, showing a substantially stronger conservation of binding energy than expected for neutral evolution [16, 27, 39].
We can now build a probabilistic model for crossspecies comparisons [16]. It is based on the joint distributions of energy pairs
P_{0}(E_{1}, E_{2}) = G_{0}(E_{2}E_{1}) P_{0}(E_{1}) (59)
under neutral evolution and
Q(E_{1}, E_{2}) = G_{ f }(E_{2}E_{1}) Q(E_{1}) (60)
under stationary selection, which are determined by the corresponding distributions in one species and the energy transition probabilities. Detailed balance of the substitution dynamics implies
i.e., the joint distributions P_{0}(E_{1}, E_{2}) and Q(E_{1}, E_{2}) must be symmetric functions of their arguments. These distributions combine into a model for pairs of aligned loci, which generalizes the singlespecies model (25) and takes the form
W(E_{1}, E_{2}) = (1  λ)P_{0}(E_{1}, E_{2}) + λQ(E_{1}, E_{2}). (62)
(This model can be extended further to include nonstationary selection.) The distribution W(E_{1}, E_{2}) with a fraction of functionality λ = 0.0018 is in excellent agreement with the count distribution W_{dat}(E_{1}, E_{2}) obtained from E. coli and S. typhimurium, as shown in fig. 8(b). The symmetry of W_{dat} is thus consistent with the underlying assumption of detailed balance. Analogous Bayesian models can be defined for more than two species related by a phylogeny. This approach has been applied to binding site prediction in bacteria [16]; a related study of several species of funghi has been reported in ref. [40].
Adaptive evolution of binding sites
What does this picture say about the adaptive evolution of transcriptional regulation in response to a newly arising selection pressure? The evolution from a genotype with marginal binding (E(a) ≈ E_{ s }) to strong binding requires only about three uphill point mutations in the fitness landscape of fig. 4(b), i.e., there is an effective fitness gain 2NΔF ≈ 3 per mutation. Hence, according to (51), the rate of uphill substitutions per locus is enhanced by a factor 2NΔF·d(a, a*) at least of order 10 over the neutral point mutation rate per nucleotide. At the same time, the downhill rate is strongly suppressed. This shows that the adaptive formation of a binding site from background sequence can indeed be a rapid mode of regulatory evolution, due to the substantial level of selection [18].
However, this mode is only efficient if adaptation can set in immediately after the selection pressure is established. In larger regulatory regions, the exact position of a binding site is often not important. We assume the initial genome contains a set of $\tilde{L}$shadow sites, i.e., positions r_{1},..., ${r}_{\tilde{L}}$ where a given sequence a would have the same regulatory effect. If one of these shadow sites has already a genotype with marginal binding, it acts as a "seed" for the onset of adaptation [41]. On the other hand, if all shadow sites of the initial genome have energy E > E_{ s }, there is typically a substantial waiting time of neutral evolution before one of them reaches the threshold energy E_{ s }. Assuming the initial genome to be entirely background sequence, it will contain at least one such seed if ${\int}_{E<{E}_{s}}{P}_{0}(E)\text{d}E}\gtrsim 1/\tilde{L$, which is a joint condition on $\tilde{L}$ and the site length ℓ: the shadow regulatory region must be long enough and binding sites must be short enough. The example shows that the evolvability of regulation imposes constraints on genome architecture [18]. Adaptive point substitution may thus be a feasible mode for the formation of a single binding site, but will hardly explain the groups of adjacent sites characteristic of eukaryotic promoters. These may originate from repeat duplication by slippage, which has recently been shown to be an efficient source of sequence innovation in intergenic regions of Drosophila.
Towards a dynamical picture of the genome
The relationship S = 2NF + const. between score and fitness is a cornerstone of the theoretical picture developed so far, which links its population genetic, bioinformatic and biophysical arches. It relates a key evolutionary variable with the statistics of genomic frequency counts. The physical binding energy is an appropriate phenotypic variable on which fitness and score depend, because molecular function is determined by binding interactions.
We have discussed this picture for transcription factor binding sites, but it can be applied more generally to functional elements in genomes. It relates the statistics of these elements in one genome with their evolutionary dynamics, which is observed in crossspecies comparisons. This dynamics is shaped by selection: The components of functional elements are coupled by a common fitness function; such fitness interactions are called epistasis. Hence, functional correlations lead to evolutionary correlations. These can be traced in the Q distribution over fixed genomes of a functional element. A more detailed statistical analysis using the statistics of polymorphisms within a population is briefly sketched below.
Thus, the picture of the genome as a system with multiple interactions has a fundamental dynamical significance. This is important since it allows us to trace functional modules from evolutionary patterns. We conclude the article with a brief outlook on functional integration of regulatory sequences at various and its dynamical implications.
Evolutionary interactions between sites
Regulatory function is often determined not by single binding sites, but jointly by a group of sites in the same regulatory region [42]. An important mechanism is binding cooperativity, i.e., the formation of a protein complex between two (or more) factors bound to their corresponding DNA sites. The binding energy of this complex has the form E = E_{1} + E_{2} + ΔE_{12}, where E_{1} and E_{2} are the energies of the factors bound individually and ΔE_{12} < 0 is the energy gain due to the proteinprotein interaction, which is of the order of a few k_{ B }T . Cooperative binding has a number of functional effects [1]:

(a)
It increases the signaltonoise ratio for the targeting of regulatory input to a specific gene, which is important in larger eukaryotic genomes, where single spurious binding sites are abundant in background sequence.

(b)
It sharpens the response of the binding probability to variations in the factor concentrations around their threshold value. This follows from the thermodynamics of two factors, which is a straightforward generalization of the case of a single factor discussed above.

(c)
It implements logical connections between regulatory input signals to a given gene. The simplest example is an AND connection between two factors, where the regulated gene is affected only if both factors are simultaneously present. This happens if the binding energies and factor concentrations are such that individual binding is weak but joint binding is strong. Larger groups of binding sites can encode a whole repertoire of more complicated logical functions [43].
Regulatory modules with several jointly acting binding sites are frequently found in eukaryotes. The functional coupling of sites in a module translates into interactions between these sites in their sequence evolution. The genomic functional element, i.e., the subset of the regulatory region on which selection acts, is the module as a whole. Its fitness F(E_{1}, E_{2}, ΔE_{12},...) is a joint function of the binding energies as the relevant phenotypic variables [18, 25]. The evolutionary dynamics under this selection allows for a large number of compensatory changes, i.e., pairs of correlated substitutions changing two binding energies such that the fitness remains constant. These lead to nucleotide frequency correlations between different sites. Such compensatory changes have indeed been observed in experiments on Drosophila promoters [44].
Siteshadow interactions
In larger regulatory regions, there is a number of shadow sites where a binding sequence a would have a similar regulatory effect as at the functional sites present. In that case, the genomic functional element contains not only the functional binding sites but also the shadow sites. Once a functional site has disappeared due to deleterious mutations, a shadow site can turn functional by adaptive evolution as described in the last section. The resulting evolutionary dynamics leads to sequence turnover with the actual binding sites present at different but functionally equivalent positions [36]. Substantial sequence turnover has been observed in a number of case studies [44–49]. Also the number of actual sites is subject to evolutionary variation since the same regulatory effect, i.e., the same fitness, can be distributed over fewer stronger or more weaker sites. With increasing number $\tilde{L}$ of shadow positions, one expects that the number of actual sites grows while individual sites get weaker [36].
Gene interactions
Evolutionary interactions are not limited to regulatory elements for the same gene. An example are gene duplications and the subsequent evolution of the daughter genes. Selection acts jointly on this pair of genes [50], which have initially identical functions, eventually leading to either loss of one of them or to subfunctionalization, which has been argued to be an important mode of genome evolution in eukaryotes [51, 52]. This process can take place by regulation, i.e., via a correlated distribution of the regulatory elements on the daughter genes. More generally, the evolution of genes in a regulatory network is correlated if their functions are coupled either in series (i.e., one gene acts on the other) or in parallel (i.e., they are part of alternative pathways for the same function). Although some regulatory networks in model organisms – e.g. the embryonic development in the sea urchin [53] – have been studied in detail, we lack a coherent view of their functional evolution to date.
Interactions and timedependent selection
The functional integration of regulation at multiple levels and the resulting fitness interactions (epistasis) imply that the selection at one genomic site is influenced by changes at other sites. A recent analysis of singlenucleotide polymorphisms and substitutions in Drosophila provides indeed evidence on a genomewide scale that selection is timedependent: at individual loci, changes in the direction of selection occur at nearly the rate of neutral evolution [54, 55]. At the same time, selection is sufficiently strong so that the adaptive response can keep up with the rate of selection changes. This rate is faster in noncoding DNA, which points towards the role of regulation in the adaptive differentiation between species. Genomic evolution emerges as a complex stochastic process, shaped jointly be the driving force of timedependent selection, fitness interactions between sites, and the ongoing background of nearneutral changes. Much more remains to be learned about the interplay of these evolutionary forces: in a large and strongly coupled system, one external signal can trigger an avalanche of subsequent compensatory responses. This dynamics seems now within reach of genomic sequence analysis.
Evolutionary innovations
Under stationary selection, functional elements are more conserved than background sequence, and the scorefitness relation quantifies the amount of conservation. But evolution is, of course, not limited to conservation. On one hand, there is typically a multitude of different genotypes yielding the same molecular function, and the evolutionary dynamics continuously plays with these alternatives. On the other hand, organisms face longterm changes of their environment, which lead to new selection pressures and a response by adaptive evolution of new functions. If regulation is to account for a large part of the diversification in higher eukaryotes, loss or gain of regulatory function should be an important mode of molecular evolution. Changes in regulatory DNA leading to new functions of gene networks have been observed [56], and it is possible to extend the statistical models described in the previous section to include evolutionary gain or loss of function of individual binding sites [16]. On a broader scale, timedependent selection and fitness couplings appear act as a major driving forces of genomic change, triggering avalanches of evolutionary innovation. Understanding this molecular basis of innovations is a major challenge for theory and experiment in the coming years. It will profoundly change our dynamical view of the genome.
References
 1.
Ptashne M, Gann A: Genes and Signals. 2002, Cold Spring Harbor Laboratory Press
 2.
Tautz D: Evolution of transcriptional regulation. Curr Opin Genet Dev. 2000, 10 (5): 575579. 10.1016/S0959437X(00)001301.
 3.
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. 10.1093/molbev/msg140. [http://dx.doi.org/10.1093/molbev/msg140]
 4.
Berg OG, Winter RB, von Hippel PH: Diffusiondriven mechanisms of protein translocation on nucleic acids. 1. Models and theory. Biochemistry. 1981, 20 (24): 69296948. 10.1021/bi00527a028.
 5.
Winter RB, von Hippel PH: Diffusiondriven mechanisms of protein translocation on nucleic acids. 2. The Escherichia coli repressoroperator interaction: equilibrium measurements. Biochemistry. 1981, 20 (24): 69486960. 10.1021/bi00527a029.
 6.
Winter RB, Berg OG, von Hippel PH: Diffusiondriven mechanisms of protein translocation on nucleic acids. 3. The Escherichia coli lac repressoroperator interaction: kinetic measurements and conclusions. Biochemistry. 1981, 20 (24): 69616977. 10.1021/bi00527a030.
 7.
von Hippel PH, Berg OG: On the specificity of DNAprotein interactions. Proc Natl Acad Sci USA. 1986, 83 (6): 16081612. 10.1073/pnas.83.6.1608.
 8.
Sarai A, Takeda Y: Lambda repressor recognizes the approximately 2fold symmetric halfoperator sequences asymmetrically. Proc Natl Acad Sci USA. 1989, 86 (17): 65136517. 10.1073/pnas.86.17.6513.
 9.
Fields DS, He Y, AlUzri AY, Stormo GD: Quantitative specificity of the Mnt repressor. J Mol Biol. 1997, 271 (2): 178194. 10.1006/jmbi.1997.1171.
 10.
Stormo GD, Fields DS: Specificity, free energy and information content in proteinDNA interactions. Trends Biochem Sci. 1998, 23 (3): 109113. 10.1016/S09680004(98)011876.
 11.
Oda M, Furukawa K, Ogata K, Sarai A, Nakamura H: Thermodynamics of specific and nonspecific DNA binding by the cMyb DNAbinding domain. J Mol Biol. 1998, 276 (3): 571590. 10.1006/jmbi.1997.1564. [http://dx.doi.org/10.1006/jmbi.1997.1564]
 12.
Omagari K, Yoshimura H, Takano M, Hao D, Ohmori M, Sarai A, Suyama A: Systematic single basepair substitution analysis of DNA binding by the cAMP receptor protein in cyanobacterium Synechocystis sp. PCC 6803. FEBS Lett. 2004, 563 (1–3): 5558. 10.1016/S00145793(04)002480. [http://dx.doi.org/10.1016/S00145793(04)002480]
 13.
Foat BC, Houshmandi SS, Olivas WM, Bussemaker HJ: Profiling conditionspecific, genomewide regulation of mRNA stability in yeast. Proc Natl Acad Sci USA. 2005, 102 (49): 1767517680. 10.1073/pnas.0503803102. [http://dx.doi.org/10.1073/pnas.0503803102]
 14.
Gerland U, Moroz JD, Hwa T: Physical constraints and functional characteristics of transcription factorDNA interaction. Proc Natl Acad Sci USA. 2002, 99 (19): 1201512020. 10.1073/pnas.192693599. [http://dx.doi.org/10.1073/pnas.192693599]
 15.
Djordjevic M, Sengupta AM, Shraiman BI: A biophysical approach to transcription factor binding site discovery. Genome Res. 2003, 13 (11): 23812390. 10.1101/gr.1271603. [http://dx.doi.org/10.1101/gr.1271603]
 16.
Mustonen V, Lässig M: Evolutionary population genetics of promoters: predicting binding sites and functional phylogenies. Proc Natl Acad Sci USA. 2005, 102 (44): 1593641. 10.1073/pnas.0505537102. [http://dx.doi.org/10.1073/pnas.0505537102]
 17.
Slutzky M, Mirny L: Kinetics of proteinDNA interaction: facilitated target location in a sequencedependent potential. Biophys J. 2004, 87 (6): 40214035. 10.1529/biophysj.104.050765.
 18.
Berg J, Willmann S, Lässig M: Adaptive evolution of transcription factor binding sites. BMC Evol Biol. 2004, 4: 4210.1186/14712148442. [http://dx.doi.org/10.1186/14712148442]
 19.
Durbin R, Eddy SR, Krogh A, Mitchison G: Biological Sequence Analysis. 1998, Cambridge University Press
 20.
Stormo GD, Hartzell GW: Identifying proteinbinding sites from unaligned DNA fragments. Proc Natl Acad Sci USA. 1989, 86 (4): 11831187. 10.1073/pnas.86.4.1183.
 21.
Hertz GZ, Stormo GD: Identifying DNA and protein patterns with statistically significant alignments of multiple sequences. Bioinformatics. 1999, 15 (7–8): 563577. 10.1093/bioinformatics/15.7.563.
 22.
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. 10.1101/gr.207502. Article published online before print in January 2002. [http://www.genome.org/cgi/content/full/12/2/298]
 23.
van Nimwegen E, Zavolan M, Rajewsky N, Siggia ED: Probabilistic clustering of sequences: inferring new bacterial regulons by comparative genomics. Proc Natl Acad Sci USA. 2002, 99 (11): 73237328. 10.1073/pnas.112690399. [http://dx.doi.org/10.1073/pnas.112690399]
 24.
Lenhard B, Sandelin A, Mendoza L, Engström P, Jareborg N, Wasserman WW: Identification of conserved regulatory elements by comparative genome analysis. J Biol. 2003, 2 (2): 1310.1186/14754924213. [http://dx.doi.org/10.1186/14754924213]
 25.
Gerland U, Hwa T: On the selection and evolution of regulatory DNA motifs. J Mol Evol. 2002, 55 (4): 386400. 10.1007/s002390022335z. [http://dx.doi.org/10.1007/s002390022335z]
 26.
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: 1910.1186/14712148319. [http://dx.doi.org/10.1186/14712148319]
 27.
Brown CT, Callan CG: Evolutionary comparisons suggest many novel cAMP response protein binding sites in Escherichia coli. Proc Natl Acad Sci USA. 2004, 101 (8): 24042409. 10.1073/pnas.0308628100.
 28.
Kimura M, Crow J: An Introduction to Population Genetics Theory. 1970, Harper & Row, New York
 29.
Kimura M: The Neutral Theory of Molecular Evolution. 1983, Cambridge University Press
 30.
Kimura M: On the probability of fixation of mutant genes in a population. Genetics. 1962, 47: 713719.
 31.
Kimura M, Ohta T: The Average Number of Generations until Fixation of a Mutant Gene in a Finite Population. Genetics. 1969, 61 (3): 763771.
 32.
Rouzine IM, Rodrigo A, Coffin JM: Transition between stochastic evolution and deterministic evolution in the presence of selection: general theory and application to virology. Microbiol Mol Biol Rev. 2001, 65: 151185. 10.1128/MMBR.65.1.151185.2001. [http://dx.doi.org/10.1128/MMBR.65.1.151185.2001]
 33.
Grün D, Lässig M: to be published.
 34.
Arndt PF, Hwa T: Identification and measurement of neighbordependent nucleotide substitution processes. Bioinformatics. 2005, 21 (10): 23222328. 10.1093/bioinformatics/bti376. [http://dx.doi.org/10.1093/bioinformatics/bti376]
 35.
Peliti L: Quasispecies evolution in general meanfield landscapes. Europhys Lett. 2002, 57: 74551. 10.1209/epl/i2002005265.
 36.
Berg J, Lässig M: Stochastic evolution of transcription factor binding sites. Biophysics (Moscow). 2003, 48 (Suppl 1): S36S44.
 37.
Halpern AL, Bruno WJ: Evolutionary distances for proteincoding sequences: modeling sitespecific residue frequencies. Mol Biol Evol. 1998, 15 (7): 910917.
 38.
Sella G, Hirsh AE: The application of statistical physics to evolutionary biology. Proc Natl Acad Sci USA. 2005, 102 (27): 95419546. 10.1073/pnas.0501865102. [http://dx.doi.org/10.1073/pnas.0501865102]
 39.
Mustonen V, Lässig M: to be published.
 40.
Moses AM, Chiang DY, Pollard DA, Iyer VN, Eisen MB: MONKEY: identifying conserved transcriptionfactor binding sites in multiple alignments using a binding sitespecific evolutionary model. Genome Biol. 2004, 5 (12): R9810.1186/gb2004512r98. [http://dx.doi.org/10.1186/gb2004512r98]
 41.
MacArthur S, Brookfield JFY: Expected rates and modes of evolution of enhancer sequences. Mol Biol Evol. 2004, 21 (6): 10641073. 10.1093/molbev/msh105. [http://dx.doi.org/10.1093/molbev/msh105]
 42.
Arnosti DN: Analysis and function of transcriptional regulatory elements: insights from Drosophila. Annu Rev Entomol. 2003, 48: 579602. 10.1146/annurev.ento.48.091801.112749. [http://dx.doi.org/10.1146/annurev.ento.48.091801.112749]
 43.
Buchler NE, Gerland U, Hwa T: On schemes of combinatorial transcription logic. Proc Natl Acad Sci USA. 2003, 100 (9): 51365141. 10.1073/pnas.0930314100. [http://dx.doi.org/10.1073/pnas.0930314100]
 44.
Ludwig MZ, Bergman C, Patel NH, Kreitman M: Evidence for stabilizing selection in a eukaryotic enhancer element. Nature. 2000, 403 (6769): 564567. 10.1038/35000615. [http://dx.doi.org/10.1038/35000615]
 45.
Ludwig MZ, Patel NH, Kreitman M: Functional analysis of eve stripe 2 enhancer evolution in Drosophila: rules governing conservation and change. Development. 1998, 125 (5): 949958.
 46.
McGregor AP, Shaw PJ, Hancock JM, Bopp D, Hediger M, Wratten NS, Dover GA: Rapid restructuring of bicoiddependent hunchback promoters within and between Dipteran species: implications for molecular coevolution. Evol Dev. 2001, 3 (6): 397407. 10.1046/j.1525142X.2001.01043.x.
 47.
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): 70314. 10.1093/molbev/msg077. [http://dx.doi.org/10.1093/molbev/msg077]
 48.
Scemama JL, Hunter M, McCallum J, Prince V, Stellwag E: Evolutionary divergence of vertebrate Hoxb2 expression patterns and transcriptional regulatory loci. J Exp Zool. 2002, 294 (3): 28599. 10.1002/jez.90009. [http://dx.doi.org/10.1002/jez.90009]
 49.
Costas J, Casares F, Vieira J: Turnover of binding sites for transcription factors involved in early Drosophila development. Gene. 2003, 310: 21520. 10.1016/S03781119(03)005560.
 50.
Wagner A: Selection and gene duplication: a view from the genome. Genome Biol. 2002, 3 (5): reviews101210.1186/gb200235reviews1012.
 51.
Lynch M, Conery JS: The evolutionary demography of duplicate genes. J Struct Funct Genomics. 2003, 3 (1–4): 3544. 10.1023/A:1022696612931.
 52.
Lynch M, Conery JS: The origins of genome complexity. Science. 2003, 302 (5649): 14014. 10.1126/science.1089370. [http://dx.doi.org/10.1126/science.1089370]
 53.
Davidson E: A view from the genome: spatial control of transcription in sea urchin development. Curr Opin Genet Dev. 1999, 9 (5): 53041. 10.1016/S0959437X(99)000131.
 54.
Mustonen V, Lässig M: Adaptations to fluctuating selection in Drosophila. Proc Natl Acad Sci USA. 2007, 104 (7): 227782. 10.1073/pnas.0607105104. [http://dx.doi.org/10.1073/pnas.0607105104]
 55.
Mustonen V, Lässig M: Sequence evolution under quenched selection fluctuations. preprint. 2006
 56.
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): e39810.1371/journal.pbio.0020398. [http://dx.doi.org/10.1371/journal.pbio.0020398]
 57.
Robison K, McGuire AM, Church GM: A comprehensive library of DNAbinding site matrices for 55 proteins applied to the complete Escherichia coli K12 genome. J Mol Biol. 1998, 284 (2): 241254. 10.1006/jmbi.1998.2160. [http://dx.doi.org/10.1006/jmbi.1998.2160]
Acknowledgements
This article has been published as part of BMC Bioinformatics Volume 8 Supplement 6, 2007: Otto Warburg International Summer School and Workshop on Networks and Regulation. The full contents of the supplement are available online at http://www.biomedcentral.com/14712105/8?issue=S6
Author information
Rights and permissions
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
About this article
Cite this article
Lässig, M. From biophysics to evolutionary genetics: statistical aspects of gene regulation. BMC Bioinformatics 8, S7 (2007). https://doi.org/10.1186/147121058S6S7
Published:
Keywords
 Binding Energy
 Effective Population Size
 Detailed Balance
 Functional Site
 Functional Locus