Volume 8 Supplement 6

## Otto Warburg International Summer School and Workshop on Networks and Regulation

# From biophysics to evolutionary genetics: statistical aspects of gene regulation

- Michael Lässig
^{1}Email author

**8(Suppl 6)**:S7

https://doi.org/10.1186/1471-2105-8-S6-S7

© Lässig; licensee BioMed Central Ltd. 2007

**Published: **27 September 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.

## Keywords

## 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, micro-RNAs 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

*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

*cis-regulatory region*of a gene, which lies just upstream of its protein-coding 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 protein-protein 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.

### Factor-DNA binding energies

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

- (a)
The single nucleotides of a binding locus

**a**= (*a*_{1},...,*a*_{ℓ}) give approximately independent contributions to the binding energy,

- (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 minimum-energy 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.

with *ε* ≈ 2*k*_{
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

*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 low-energy 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 one-dimensional diffusion along the DNA backbone and three-dimensional diffusion in the surrounding medium. This proves more efficient than purely one- or three-dimensional 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 low-energy 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

*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

*r*=

*r*

_{ f }, must have a low specific binding energy

*E*=

*E*(

*r*

_{ f }). We now single out this position and write

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

*r*can be empty or be occupied either by an unspecifically or by a specifically bound factor. Using a chemical potential

*σ*, the many-factor partition function can hence be written as

*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 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 single-molecule 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 signal-to-noise ratio of regulatory interactions,

*F*_{0} ≳ *E**, (11)

- (a)
In a random genome, there must be at most a number of order one minimum-energy 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)

*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 minimum-energy site, there are ℓ suboptimal sites of Hamming distance 1 from the minimum-energy sequence. These must not suppress the binding to the minimum-energy site, i.e.,

exp(-

*E**/*k*_{ B }*T*) ≳ ℓ exp [-(*E** +*ε*)/*k*_{ B }*T*] (13)

*ε*/

*k*

_{ B }

*T*≳ log ℓ ≈ 2 - 3.

- (c)
Finally, the unspecific binding in the entire genome must not suppress the specific binding to a minimum-energy 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 single-molecule 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 minimum-energy 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 single-molecule 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 signal-to-noise 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

*p*

_{0}(

*a*) (

*a*=

*A*,

*C*,

*G*,

*T*). Hence, the probability of finding a given sequence has the factorized form

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

*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 cross-species 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(**a**|*m*) prob(*m*) = prob(*m*|**a**) prob(**a**) (23)

**a**) = ∑

_{ m }prob(

**a**|

*m*)prob(

*m*). We can solve for the conditional probability of the model for given data

**a**,

*ρ*

_{ f }(

**a**) = prob(

*m*= 1|

**a**), 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*).

*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

*s*of non-overlapping 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 no-overlap constraint

*r*

_{ν+1}≥

*r*

_{ ν }+ ℓ for

*ν*= 1,...,

*s*- 1). The full sequence distribution in a segment of length

*L*has the form

*Z*is a normalization factor, $\tilde{\lambda}$ =

*λ*+

*O*(

*λ*

^{2}) is a weight factor for each functional locus (the negligible correction terms originate from the no-overlap 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*_{r-1}(*a*_{1},..., *a*_{r-1}) + $\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.

**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 no-overlap constraint is then ${r}_{f}^{\ast}$ = {*r|b*(*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 non-overlapping 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 self-contained 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.

*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 time-dependent 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).

*F*

_{ ab }to be constant over the time of observation, the solution of eq. (32) is the evolutionary trajectory

*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

*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,

*χ*

_{ 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 large-*N* 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 so-called *effective population size*, which can be inferred from genome data and is in general smaller than the actual population size.

*x*,

*χ*

_{ x }(

*t*) = (

*∂x*/

*∂N*

_{ a })

*χ*

_{ a }(

*t*) + (

*∂x*/

*∂N*

_{ b })

*χ*

_{ b }(

*t*) are Gaussian random variables with zero mean and

*x*(

*t*) as shown in fig. 6(b). To capture their statistics, we convert the Langevin equation (36) into a Fokker-Planck equation for the probability distribution of the genotype composition [28, 30],

*x*-dependent diffusion "constant"

*x*(1 -

*x*)/2

*N*, 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 Fokker-Planck equation (38) describes diffusion in the interval (0, 1) with

*absorbing boundaries*. There is a family of stationary states

*x*) = (1 -

*φ*)

*δ*(

*x*) +

*φδ*(1 -

*x*),

*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

*F*

_{ ab }, the stationary solution

*φ*(

*x*

_{0}) = lim

_{t→∞}$\mathcal{\text{P}}$(

*x*= 1,

*t*|

*x*

_{0},

*t*

_{0}) has the form [28, 30]

which for near-neutral 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

**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,

*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*,

*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}(N-1))$[33].

- (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 so-called*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)

*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
}+ ...).

*μ*

_{a→b}are referred to as neutral mutation rates. For strong selection (

*N*|Δ

*F*

_{ ab }| ≫ 1 ≫ |Δ

*F*

_{ ab }|), eq. (49) takes the asymptotic forms

*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

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

*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 score-fitness 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 [2*NF* (**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 [2*NF* + *S*_{0}(*F*) + const.] (57)

For binding site sequences, we have a similar projection on the binding energy, *Q*(*E*) = exp [2*NF*(*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 2*NF* - *S*_{0} is dominated by the fitness term, and *Q*(*F*) takes appreciable values only at points of near-maximal fitness, i.e., where *F*_{
max
}- *F* ≲ 1/2*N*. 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 2*N* 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 log-likelihood 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**) = 2*NF*(**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 2*NF*(*E*) for CRP binding sites is thus given by fig. 4(b)[16]. The fitness is constant in the no-binding 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 2*NF* 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

*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 2

*NF*(

*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 cross-species 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)

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 single-species 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 non-stationary 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 2*NΔF* ≈ 3 per mutation. Hence, according to (51), the rate of uphill substitutions per locus is enhanced by a factor 2*NΔ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* = 2*NF* + 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 cross-species 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

*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 protein-protein 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 signal-to-noise 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].

### Site-shadow 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 time-dependent 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 single-nucleotide polymorphisms and substitutions in *Drosophila* provides indeed evidence on a genome-wide scale that selection is time-dependent: 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 non-coding 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 time-dependent selection, fitness interactions between sites, and the ongoing background of near-neutral 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 score-fitness 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 long-term 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, time-dependent 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.

## Declarations

### 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/1471-2105/8?issue=S6

## Authors’ Affiliations

## References

- Ptashne M, Gann A: Genes and Signals. 2002, Cold Spring Harbor Laboratory PressGoogle Scholar
- Tautz D: Evolution of transcriptional regulation. Curr Opin Genet Dev. 2000, 10 (5): 575-579. 10.1016/S0959-437X(00)00130-1.View ArticlePubMedGoogle Scholar
- 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): 1377-1419. 10.1093/molbev/msg140. [http://dx.doi.org/10.1093/molbev/msg140]View ArticlePubMedGoogle Scholar
- Berg OG, Winter RB, von Hippel PH: Diffusion-driven mechanisms of protein translocation on nucleic acids. 1. Models and theory. Biochemistry. 1981, 20 (24): 6929-6948. 10.1021/bi00527a028.View ArticlePubMedGoogle Scholar
- Winter RB, von Hippel PH: Diffusion-driven mechanisms of protein translocation on nucleic acids. 2. The Escherichia coli repressor-operator interaction: equilibrium measurements. Biochemistry. 1981, 20 (24): 6948-6960. 10.1021/bi00527a029.View ArticlePubMedGoogle Scholar
- Winter RB, Berg OG, von Hippel PH: Diffusion-driven mechanisms of protein translocation on nucleic acids. 3. The Escherichia coli lac repressor-operator interaction: kinetic measurements and conclusions. Biochemistry. 1981, 20 (24): 6961-6977. 10.1021/bi00527a030.View ArticlePubMedGoogle Scholar
- von Hippel PH, Berg OG: On the specificity of DNA-protein interactions. Proc Natl Acad Sci USA. 1986, 83 (6): 1608-1612. 10.1073/pnas.83.6.1608.PubMed CentralView ArticlePubMedGoogle Scholar
- Sarai A, Takeda Y: Lambda repressor recognizes the approximately 2-fold symmetric half-operator sequences asymmetrically. Proc Natl Acad Sci USA. 1989, 86 (17): 6513-6517. 10.1073/pnas.86.17.6513.PubMed CentralView ArticlePubMedGoogle Scholar
- Fields DS, He Y, Al-Uzri AY, Stormo GD: Quantitative specificity of the Mnt repressor. J Mol Biol. 1997, 271 (2): 178-194. 10.1006/jmbi.1997.1171.View ArticlePubMedGoogle Scholar
- Stormo GD, Fields DS: Specificity, free energy and information content in protein-DNA interactions. Trends Biochem Sci. 1998, 23 (3): 109-113. 10.1016/S0968-0004(98)01187-6.View ArticlePubMedGoogle Scholar
- Oda M, Furukawa K, Ogata K, Sarai A, Nakamura H: Thermodynamics of specific and non-specific DNA binding by the c-Myb DNA-binding domain. J Mol Biol. 1998, 276 (3): 571-590. 10.1006/jmbi.1997.1564. [http://dx.doi.org/10.1006/jmbi.1997.1564]View ArticlePubMedGoogle Scholar
- Omagari K, Yoshimura H, Takano M, Hao D, Ohmori M, Sarai A, Suyama A: Systematic single base-pair substitution analysis of DNA binding by the cAMP receptor protein in cyanobacterium Synechocystis sp. PCC 6803. FEBS Lett. 2004, 563 (1–3): 55-58. 10.1016/S0014-5793(04)00248-0. [http://dx.doi.org/10.1016/S0014-5793(04)00248-0]View ArticlePubMedGoogle Scholar
- Foat BC, Houshmandi SS, Olivas WM, Bussemaker HJ: Profiling condition-specific, genome-wide regulation of mRNA stability in yeast. Proc Natl Acad Sci USA. 2005, 102 (49): 17675-17680. 10.1073/pnas.0503803102. [http://dx.doi.org/10.1073/pnas.0503803102]PubMed CentralView ArticlePubMedGoogle Scholar
- Gerland U, Moroz JD, Hwa T: Physical constraints and functional characteristics of transcription factor-DNA interaction. Proc Natl Acad Sci USA. 2002, 99 (19): 12015-12020. 10.1073/pnas.192693599. [http://dx.doi.org/10.1073/pnas.192693599]PubMed CentralView ArticlePubMedGoogle Scholar
- Djordjevic M, Sengupta AM, Shraiman BI: A biophysical approach to transcription factor binding site discovery. Genome Res. 2003, 13 (11): 2381-2390. 10.1101/gr.1271603. [http://dx.doi.org/10.1101/gr.1271603]PubMed CentralView ArticlePubMedGoogle Scholar
- Mustonen V, Lässig M: Evolutionary population genetics of promoters: predicting binding sites and functional phylogenies. Proc Natl Acad Sci USA. 2005, 102 (44): 15936-41. 10.1073/pnas.0505537102. [http://dx.doi.org/10.1073/pnas.0505537102]PubMed CentralView ArticlePubMedGoogle Scholar
- Slutzky M, Mirny L: Kinetics of protein-DNA interaction: facilitated target location in a sequence-dependent potential. Biophys J. 2004, 87 (6): 4021-4035. 10.1529/biophysj.104.050765.View ArticleGoogle Scholar
- Berg J, Willmann S, Lässig M: Adaptive evolution of transcription factor binding sites. BMC Evol Biol. 2004, 4: 42-10.1186/1471-2148-4-42. [http://dx.doi.org/10.1186/1471-2148-4-42]PubMed CentralView ArticlePubMedGoogle Scholar
- Durbin R, Eddy SR, Krogh A, Mitchison G: Biological Sequence Analysis. 1998, Cambridge University PressView ArticleGoogle Scholar
- Stormo GD, Hartzell GW: Identifying protein-binding sites from unaligned DNA fragments. Proc Natl Acad Sci USA. 1989, 86 (4): 1183-1187. 10.1073/pnas.86.4.1183.PubMed CentralView ArticlePubMedGoogle Scholar
- Hertz GZ, Stormo GD: Identifying DNA and protein patterns with statistically significant alignments of multiple sequences. Bioinformatics. 1999, 15 (7–8): 563-577. 10.1093/bioinformatics/15.7.563.View ArticlePubMedGoogle Scholar
- Rajewsky N, Socci ND, Zapotocky M, Siggia ED: The evolution of DNA regulatory regions for proteo-gamma bacteria by interspecies comparisons. Genome Res. 2002, 12 (2): 298-308. 10.1101/gr.207502. Article published online before print in January 2002. [http://www.genome.org/cgi/content/full/12/2/298]PubMed CentralView ArticlePubMedGoogle Scholar
- 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): 7323-7328. 10.1073/pnas.112690399. [http://dx.doi.org/10.1073/pnas.112690399]PubMed CentralView ArticlePubMedGoogle Scholar
- 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): 13-10.1186/1475-4924-2-13. [http://dx.doi.org/10.1186/1475-4924-2-13]PubMed CentralView ArticlePubMedGoogle Scholar
- Gerland U, Hwa T: On the selection and evolution of regulatory DNA motifs. J Mol Evol. 2002, 55 (4): 386-400. 10.1007/s00239-002-2335-z. [http://dx.doi.org/10.1007/s00239-002-2335-z]View ArticlePubMedGoogle Scholar
- 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-10.1186/1471-2148-3-19. [http://dx.doi.org/10.1186/1471-2148-3-19]PubMed CentralView ArticlePubMedGoogle Scholar
- 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): 2404-2409. 10.1073/pnas.0308628100.PubMed CentralView ArticlePubMedGoogle Scholar
- Kimura M, Crow J: An Introduction to Population Genetics Theory. 1970, Harper & Row, New YorkGoogle Scholar
- Kimura M: The Neutral Theory of Molecular Evolution. 1983, Cambridge University PressView ArticleGoogle Scholar
- Kimura M: On the probability of fixation of mutant genes in a population. Genetics. 1962, 47: 713-719.PubMed CentralPubMedGoogle Scholar
- Kimura M, Ohta T: The Average Number of Generations until Fixation of a Mutant Gene in a Finite Population. Genetics. 1969, 61 (3): 763-771.PubMed CentralPubMedGoogle Scholar
- 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: 151-185. 10.1128/MMBR.65.1.151-185.2001. [http://dx.doi.org/10.1128/MMBR.65.1.151-185.2001]PubMed CentralView ArticlePubMedGoogle Scholar
- Grün D, Lässig M: to be published.Google Scholar
- Arndt PF, Hwa T: Identification and measurement of neighbor-dependent nucleotide substitution processes. Bioinformatics. 2005, 21 (10): 2322-2328. 10.1093/bioinformatics/bti376. [http://dx.doi.org/10.1093/bioinformatics/bti376]View ArticlePubMedGoogle Scholar
- Peliti L: Quasispecies evolution in general mean-field landscapes. Europhys Lett. 2002, 57: 745-51. 10.1209/epl/i2002-00526-5.View ArticleGoogle Scholar
- Berg J, Lässig M: Stochastic evolution of transcription factor binding sites. Biophysics (Moscow). 2003, 48 (Suppl 1): S36-S44.Google Scholar
- Halpern AL, Bruno WJ: Evolutionary distances for protein-coding sequences: modeling site-specific residue frequencies. Mol Biol Evol. 1998, 15 (7): 910-917.View ArticlePubMedGoogle Scholar
- Sella G, Hirsh AE: The application of statistical physics to evolutionary biology. Proc Natl Acad Sci USA. 2005, 102 (27): 9541-9546. 10.1073/pnas.0501865102. [http://dx.doi.org/10.1073/pnas.0501865102]PubMed CentralView ArticlePubMedGoogle Scholar
- Mustonen V, Lässig M: to be published.Google Scholar
- Moses AM, Chiang DY, Pollard DA, Iyer VN, Eisen MB: MONKEY: identifying conserved transcription-factor binding sites in multiple alignments using a binding site-specific evolutionary model. Genome Biol. 2004, 5 (12): R98-10.1186/gb-2004-5-12-r98. [http://dx.doi.org/10.1186/gb-2004-5-12-r98]PubMed CentralView ArticlePubMedGoogle Scholar
- MacArthur S, Brookfield JFY: Expected rates and modes of evolution of enhancer sequences. Mol Biol Evol. 2004, 21 (6): 1064-1073. 10.1093/molbev/msh105. [http://dx.doi.org/10.1093/molbev/msh105]View ArticlePubMedGoogle Scholar
- Arnosti DN: Analysis and function of transcriptional regulatory elements: insights from Drosophila. Annu Rev Entomol. 2003, 48: 579-602. 10.1146/annurev.ento.48.091801.112749. [http://dx.doi.org/10.1146/annurev.ento.48.091801.112749]View ArticlePubMedGoogle Scholar
- Buchler NE, Gerland U, Hwa T: On schemes of combinatorial transcription logic. Proc Natl Acad Sci USA. 2003, 100 (9): 5136-5141. 10.1073/pnas.0930314100. [http://dx.doi.org/10.1073/pnas.0930314100]PubMed CentralView ArticlePubMedGoogle Scholar
- Ludwig MZ, Bergman C, Patel NH, Kreitman M: Evidence for stabilizing selection in a eukaryotic enhancer element. Nature. 2000, 403 (6769): 564-567. 10.1038/35000615. [http://dx.doi.org/10.1038/35000615]View ArticlePubMedGoogle Scholar
- 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): 949-958.PubMedGoogle Scholar
- McGregor AP, Shaw PJ, Hancock JM, Bopp D, Hediger M, Wratten NS, Dover GA: Rapid restructuring of bicoid-dependent hunchback promoters within and between Dipteran species: implications for molecular coevolution. Evol Dev. 2001, 3 (6): 397-407. 10.1046/j.1525-142X.2001.01043.x.View ArticlePubMedGoogle Scholar
- 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): 703-14. 10.1093/molbev/msg077. [http://dx.doi.org/10.1093/molbev/msg077]View ArticlePubMedGoogle Scholar
- 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): 285-99. 10.1002/jez.90009. [http://dx.doi.org/10.1002/jez.90009]View ArticlePubMedGoogle Scholar
- Costas J, Casares F, Vieira J: Turnover of binding sites for transcription factors involved in early Drosophila development. Gene. 2003, 310: 215-20. 10.1016/S0378-1119(03)00556-0.View ArticlePubMedGoogle Scholar
- Wagner A: Selection and gene duplication: a view from the genome. Genome Biol. 2002, 3 (5): reviews1012-10.1186/gb-2002-3-5-reviews1012.PubMed CentralView ArticlePubMedGoogle Scholar
- Lynch M, Conery JS: The evolutionary demography of duplicate genes. J Struct Funct Genomics. 2003, 3 (1–4): 35-44. 10.1023/A:1022696612931.View ArticlePubMedGoogle Scholar
- Lynch M, Conery JS: The origins of genome complexity. Science. 2003, 302 (5649): 1401-4. 10.1126/science.1089370. [http://dx.doi.org/10.1126/science.1089370]View ArticlePubMedGoogle Scholar
- Davidson E: A view from the genome: spatial control of transcription in sea urchin development. Curr Opin Genet Dev. 1999, 9 (5): 530-41. 10.1016/S0959-437X(99)00013-1.View ArticlePubMedGoogle Scholar
- Mustonen V, Lässig M: Adaptations to fluctuating selection in Drosophila. Proc Natl Acad Sci USA. 2007, 104 (7): 2277-82. 10.1073/pnas.0607105104. [http://dx.doi.org/10.1073/pnas.0607105104]PubMed CentralView ArticlePubMedGoogle Scholar
- Mustonen V, Lässig M: Sequence evolution under quenched selection fluctuations. preprint. 2006Google Scholar
- Gasch AP, Moses AM, Chiang DY, Fraser HB, Berardini M, Eisen MB: Conservation and evolution of cis-regulatory systems in ascomycete fungi. PLoS Biol. 2004, 2 (12): e398-10.1371/journal.pbio.0020398. [http://dx.doi.org/10.1371/journal.pbio.0020398]PubMed CentralView ArticlePubMedGoogle Scholar
- Robison K, McGuire AM, Church GM: A comprehensive library of DNA-binding site matrices for 55 proteins applied to the complete Escherichia coli K-12 genome. J Mol Biol. 1998, 284 (2): 241-254. 10.1006/jmbi.1998.2160. [http://dx.doi.org/10.1006/jmbi.1998.2160]View ArticlePubMedGoogle Scholar

## Copyright

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.