{'fraction_non_alphanumeric': 0.0392391870766024, 'fraction_numerical': 0.014156678825777315, 'mean_word_length': 4.32079482439926, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the Selection and Evolution of Regulatory DNA Motifs

The mutation and selection of regulatory DNA sequences is presented as an ideal model system of molecular evolution where genotype, phenotype, and fitness can be explicitly and independently characterized. In this theoretical study, we construct an explicit model for the evolution of regulatory sequences, making use of the known biophysics of the binding of regulatory proteins to DNA sequences, under the assumption that fitness of a sequence depends only on its binding affinity to the regulatory protein. The model is confined to the mean field (i.e., infinite population size) limit. Using realistic values for all parameters, we determine the minimum fitness advantage needed to maintain a binding sequence, demonstrating explicitly the ``error threshold'' below which a binding sequence cannot survive the accumulated effect of mutation over long time. The commonly observed ``fuzziness'' in binding motifs arises naturally as a consequence of the balance between selection and mutation in our model. In addition, we devise a simple model for the evolution of multiple binding sequences in a given regulatory region. We find the number of evolutionarily stable binding sequences to increase in a step-like fashion with increasing fitness advantage, if multiple regulatory proteins can synergistically enhance gene transcription. We discuss possible experimental approaches to resolve open questions raised by our study.


I. INTRODUCTION
The regulation of gene expression involves many different proteins known as transcription factors which bind passively to specific sites on the genomic DNA (see, e.g. Gerhart and Kirschner (1997)). In bacteria, each such site (called 'operator') typically consists of a contiguous sequence of 20 − 30 nucleotides which binds a specific transcription factor with much higher affinity than would a random DNA sequence of comparable length (Stormo and Fields, 1998;von Hippel, 1979). Known examples of different operators for the same factor usually differ from the maximum affinity binding sequence in a number of positions, typically in as many as 20% to 30% of the significant positions that contribute most to the specificity of the interaction. The ensemble of viable binding sequences is collectively referred to as the binding "motif" for a factor; its "fuzziness" creates a difficult computational problem for the prediction of binding sites via informatic methods (see, e.g., Lawrence et al. (1993); Stormo and Hartzell (1989); and references therein). In many known cases, a single regulatory region contains multiple operators for the same factor, each of which deviates from the maximum affinity binding sequence.
Why are the motifs fuzzy? One possible scenario is that the binding affinity of each operator is tuned evolutionarily to maximize the function of each regulated gene or operon. An alternative scenario is that the func-tion is insensitive to the detail of the binding affinity as long as it is above some threshold. In the former case, fuzziness in the binding arises due to the particular distribution of functional requirment. In the latter case, binding sequences in different regulatory regions are deemed "equal", and fuzziness results from maximizing the sequence "entropy". While anecdotal examples of both cases are known, understanding whether either case dominates in biology is of interest not only for its own sake, but also very important for the choice of proper informatics tools for motif finding. In this paper, we describe a detailed theoretical study of the latter case from an evolutionary perspective, recognizing that as with any other portion of the genome, the binding sequences are subject to the opposing forces of mutation and selection over evolutionary timescales. In particular, we address the quantitative question of how large a selective advantage the presence of a binding motif needs to provide, in order to guarantee its survival against mutations, and how large an advantage before multiple motifs are justified. To make the study concrete and explicit, we will mostly confine the discussion to gene regulation in bacteria or phages, and focus on the binding of one specific transcription factor to its operator(s) in the regulatory region of one specific gene or operon. We will not treat the interactions among different factors, since in bacteria such as E. coli, the majority of genes are regulated by a single factor (Gralla and Collado-Vides, 1996).
Another motivation for our study is that the evolution of transcription factor binding motifs seems to be a well-suited starting point for an attempt to establish a link between the microscopic molecular mechanisms in the cell and the "macroscopic" principles of evolution: In general, the most important ingredient in an evolutionary study is to relate the genotype on which mutation acts to the fitness of the organism on which selection acts through some quantifiable phenotype. This relation is particularly simple for the operator binding problem at hand, where a natural choice of the phenotype is the binding probability of the transcription factor to the operator. Regardless of whether the factor acts as an activator by attracting a polymerase to transcribe the gene, or as a repressor to block transcription, it can function only when it is bound to its operator. The fraction of time an operator is occupied in equilibrium is given by the binding probability. To regulate the transcription of the gene, e.g., in reaction to a change in the environmental conditions or in order to trigger a different phase of the cell cycle, the cell changes the factor-operator binding probability by varying the concentration of the (activated) factor inside the cell. The concentration may vary from practically zero in the "OFF-state" to typically several hundred copies per cell in the "ON-state". We will make the reasonable (but critical) assumption that the fitness gain an operator contributes to the organism depends solely on the binding probability P in the ON-state, with the value of P itself determined by the actual sequence of the operator through the binding energetics.
For a few exemplary transcription factors, the variation in binding affinity upon mutation of the binding sequence has been studied in great detail experimentally (Fields et al., 1997;Oda et al., 1998;Sarai and Takeda, 1989;Takeda et al., 1989). In particular, Fields, Stormo, and coworkers have shown for the case of the mnt repressor that its binding (free) energy is approximately a sum over independent contributions from each of the nucleotide positions in the binding sequence (Fields et al., 1997). Typically, only 10 ∼ 15 positions in a binding site have a strong preference for specific nucleotides, while the other positions do not contribute significantly to the binding energy. Known binding sequences display a fuzziness of up to 3 ∼ 4 mismatches in these significant positions. A useful simplified 'twostate model' for transcription factor binding is obtained by taking only the significant bases into account and assigning to each of them the same binding energy ε, i.e. a match (to the optimal binding sequence) is favored by an energy difference ε over a mismatch. This model, introduced long ago by von Hippel and Berg (1986), takes into account the effect of sequence-specific binding by a single parameter ε, without reference to detailed binding energies which have not yet been measured for most transcription factors.
Based on the two-state model and our assumption on the contribution of the binding of the factor towards fitness, we construct an explicit theory for the evolution of the binding sequences. Within the mean-field approach originally proposed by Eigen in the context of quasispecies evolution (Eigen, 1971;Eigen et al., 1989), we characterize the balance between the opposing forces of selection and mutation quantitatively. We determine the critical selection pressure needed to keep a motif from mutating away, and show how the fuzziness in the motifs arises naturally above the selection threshold. We further apply the theory to investigate the frequently observed occurrence of multiple motifs in a given regulatory region, and elaborate on various plausible causes. Towards the end, we will provide extended discussions on experimental approaches to pursue the open questions suggested by this study.

II. MODEL AND EQUATIONS
We focus on the operator sequence located in the regulatory region of a gene of interest. By assumption, this gene is regulated by a single transcription factor. Let S = {S 1 , S 2 , ..., S L } denote the L significant nucleotides of the operator which specify transcription factor binding. We will keep the alphabet size, A, as a variable in our equations, since, as we will see below, this facilitates the intuitive understanding of certain dependencies; however A = 4 and S i ∈ {A, C, G, T} is the only case of interest here. To describe the evolution of S in a population of bacteria or phages, we need to specify the action of selection and mutation.

Selection mechanism
It should be clear that gene regulation is only needed in the presence of a changing cellular state, triggered either internally, e.g. cell cycle, or externally through a change in the environment. Hence to study the fitness of a regulatory mechanism, we must invoke at least two different states. Selection arises when the growth rate of an organism depends on the probability P S that the factor binds to the sequence S in the state that prefers factor binding (the "ON-state"). For the sake of concreteness, let us consider an environment that oscillates between two states. We assume that in State A (the ONstate), the environment induces a certain concentration of activated factors, say on average N T F factors per cell (these may either be produced upon entering State A, or preexisting factors may be activated for binding by inducers that cause an allosteric transition). Let the growth rate or "fitness" of the organism in this state be φ A if the factor is never bound (binding motif not present), and φ A + δφ A > φ A if the factor is always bound (see Table 1). When the environment is in State B (the OFFstate), let the fitness be φ B if the factor is never bound and φ B − δφ B < φ B if the factor is always bound. In the following we will assume that the concentration of activated factors in State B is practically zero, so that the operator is never occupied (hence the parameter δφ B does not enter our model).
An example for the general situation described above could be the binding of the lac-repressor to its operator State A State B factor unbound φA φB factor bound φA + δφA φB − δφB in the lac-operon of E. coli. In this case, the ON-state would be the glucose-rich environment, and the OFFstate would be the glucose-poor and lactose-rich environment. φ A,B would be the growth rate of E. coli in the two environments in the absence of the lac-repressor. δφ A would be the increment in fitness when the unnecessary lac-operon is turned off, and φ B − δφ B ≈ 0 is the deleterious situation when lactose is present as the main source of sugar but the lac-operon is not operative due to the undesirable binding of the repressor.
In this study, we will discuss mainly the time-averaged effect over evolutionary time scales, which are much larger than the time scales of cellular or environmental fluctuations. We choose τ / ln 2 as our unit of time, with τ denoting the average generation time in the absence of the factor, so that the time averaged growth rate there can be set to 1. We assume that the cell can quickly adjust the cellular concentration of the factor 1 so that the operator with sequence S is occupied with probability P S in the ON-state and never occupied in the OFF-state. It is then plausible to assume that the time averaged growth rate Φ S depends linearly on P S (see also the discussion in the section entitled 'Selection threshold and fuzzy motifs'), Here, α is a dimensionless parameter which characterizes the selection pressure on the binding sequence S. In the limit α ≪ 1, there is hardly any selection pressure on the sequence at all; the opposite limit α → ∞ corresponds to the case where the failure of the factor-operator binding is lethal to the organism. If the fraction of time the bacteria population encounters environment A is f , the selection pressure becomes 1 In the present article, we do not consider the 'search problem' of how a transcription factor locates its operator among millions of other sites on the DNA (see Berg et al. (1981); ;  for a thorough experimental and theoretical investigation of this problem, and Gerland, Moroz, and Hwa (2001, submitted for publication) for a bound on the protein-DNA interaction parameters that results from the requirement of reasonable search times). Rather we treat protein-DNA binding as an equilibrium process characterized only by the binding probability. This is justified by the fact that the search time is typically on the order of 1 min, which is much smaller than the characteristic time scale of changes in gene expression.
In an experiment, α can be adjusted according to Eq. (2) by changing the fraction of time f the ON-state is presented. Below, we will investigate the statistics of the selected sequence S for a wide range of α's.

Mutation process
We consider only single-nucleotide substitutions, and focus on mutations in the binding sequence S, assuming that the net result of mutation and selection on the rest of the genome gives the overall background fitness of 1 (with our time unit of τ / ln 2). Furthermore, we neglect the difference between transversions and transitions, and assume a constant rate ν 0 at which a base mutates into any other base. The total mutation rate of a site of length L is then ν = ν 0 L. For bacteria such as E. Coli, ν 0 is on the order of 10 −9 under normal conditions and hence ν ∼ 10 −8 . The mutation rate is larger for RNA viruses which rely on the less accurate reverse transcriptase for replication. For that case, ν 0 is in the range 10 −5 to 10 −4 and hence ν ∼ 10 −4 − 10 −3 .

Binding probability
As mentioned above, the binding (free) energy E S of the transcription factor to the binding sequence S is given, to a good approximation, by a sum over independent contributions from each nucleotide at the L significant positions (Fields et al., 1997;Stormo and Fields, 1998), i.e. E S = L i=1 ε i (S i ). Each of these positions typically prefers a particular nucleotide by a binding energy of several k B T 's (we exclude from S those positions which contribute only a fraction of k B T towards the total binding energy). Furthermore we adopt the 'two-state model' (Berg and von Hippel, 1987;von Hippel and Berg, 1986), by assuming that each ε i (S i ) can only take on two possible values, 0 if S i matches the preferred base S * i , or ε > 0 for a mismatch, i.e. ε i (S i ) = ε · (1 − δ Si,S * i ). The binding energy of a site S is then only a function of the number of mismatches, or Hamming distance r S = | S − S * |, from the optimal binding sequence S * , i.e. 2 , Given its binding energy, the average occupancy of a site 2 Note that the approximate linear relationship between the binding energy and the number of mismatches, Eq. (3), breaks down when E reaches a certain non-specific binding energy Ens . Beyond this point, the binding energy remains constant. However, the expression (4) nevertheless provides an useful description of the binding probability over the entire range 0 ≤ r ≤ L, since P (r) is essentially zero when E = Ens.
is determined by equilibrium thermodynamics. Since a binding site can only be occupied or unoccupied (but not multiply occupied), its binding probability P S = P (r S ) is given by a Fermi function, which is also known as Arrhenius function, see, e.g., Atkins (1998). Here, µ = ε r 0 is the chemical potential for the transcription factors in the ON-state (this function is plotted in Fig. 1 (right) with realistic parameter values). Note that r 0 corresponds to the number of mismatches for which the probability of binding is 50%.
In total, we are left with three dimensionless parameters for the two-state model of protein-DNA binding: L, ε/k B T , and r 0 . As mentioned before, the number of significant positions in a binding site is typically in the range 10 ≤ L ≤ 15. By inspection of the known binding energies for exemplary transcription factors (Fields et al., 1997;Oda et al., 1998;Sarai and Takeda, 1989;Takeda et al., 1989), we find the mean specificity of the significant sites to be typically ε = 1 ∼ 3 k B T . In (Gerland, Moroz, and Hwa, 2001, submitted for publication), it is argued on rather general ground that this is actually the optimal range of ε for the transcription factors. The chemical potential µ depends directly on the average number of factors N TF in the cell; the work of Gerland, Moroz, and Hwa (2001, submitted for publication) suggests that it can be approximated by µ ≈ µ 0 + k B T ln[N TF ], where µ 0 represents the binding free energy of a single factor to the rest of the genome. For those factors whose binding energies ε i (S i ) have been measured, we find µ 0 ≈ −0.8k B T (mnt) and µ 0 ≈ −1.9k B T (lambda-repressor and cro). Hence, µ ≈ k B T ln[N TF ]; see Fig. 1 for details. For ε = 2k B T and N TF = 50 ∼ 5000, we get r 0 = µ/ε = 2 ∼ 4.3. Clearly, r 0 is the parameter that we have the least information about; but we see that it has a limited range, and in any case, most of our qualitative conclusions will be insensitive to the specific choice of r 0 . [Note that the above analysis is for factors that have a binding site only in a single regulatory region. For those factors which are global regulators and have many operators located throughout the genome (e.g., the factor CRP in E. coli), the number N TF above needs to be appropriately adjusted by the number of operators (Sengupta et al., 2002).]

Evolution equation
In this study, we will focus on the steady-state properties of the mutation/selection process defined above. For large population size and close to the steady-state, we may consider only the dynamics of the average population and neglect fluctuations due to the discreteness of the individual organisms. We denote the average number of individuals at time t with binding sequence S by  (Fields et al., 1997) on the Salmonella genome, the cro energy matrix , and the lambda repressor (cI) energy matrix  on the E. coli genome. Here, we assumed that all of the factors in the cell are bound to DNA. The dashed line delineates the slope kBT which would be expected from the considerations of Gerland, Moroz, and Hwa (2001, submitted). Right: binding probability versus the number of mismatches from the best binding sequence, according to Eq. (4) with realistic parameters.
The first term on the right hand side describes the mutational flow into N S from all sequences S ′ that are a single-nucleotide mutation away, while the second describes the reverse process. The third term represents the (time-averaged) selection/amplification process. Eq. (5) is similar to the "para-muse model" considered in a different context by Baake et al. (1997). Since the fitness function (1) depends on the sequence S only through the binding probability P S , which depends only on the number of mismatches r according to Eq. (4), it is advantageous to introduce a 'radial distribution' N (r, t) in the (discrete) Hamming distance space (Nowak and Schuster, 1989), With φ(r S ) ≡ Φ S denoting the 'radial fitness' function, the evolution equation for N (r, t) becomes where is a mesa-shaped fitness landscape. Eq. (7) is obtained by observing that there are (L − r)(A − 1) ways to mutate a site with r mismatches into a site with r+1 mismatches, r ways to mutate it into a site with r − 1 mismatches, and r(A − 2) ways to mutate a site without changing the number of mismatches. We will characterize the predictions of our model by numerically integrating the discrete radial evolution equation (7) using the set of realistic parameters given above. However, to gain insight about the qualitative behavior of the model, we also analyze the continuumspace evolution equation obtained in the limit of large L, where we used n(r, t) and ϕ(r) to denote the continuum generalization of the functions N (r, t) and φ(r) respectively. Note that the mutational dynamics is locally conservative, with a local current j(r, t) = D(r) ∂ r n(r, t) − v(r) n(r, t). The appropriate boundary conditions are j(0, t) = 0 and j(L, t) = 0. The continuous radial evolution equation (9) reduces the evolutionary dynamics to a simple one-dimensional drift-diffusion equation, where the 'diffusion coefficient' D(r) and the 'drift velocity' v(r) are explicitly given by Note that both D and v are proportional to the overall mutation rate ν = ν 0 L, with only weak dependence on r for r ≪ L. The drift velocity drives the distribution away from the optimal binding site at r = 0, simply because the number of sequences with a fixed number of mismatches r increases quickly with r. This purely entropic bias changes sign at (A − 1)L/A, which is the average number of mismatches in a random binding sequence. Also note that the r-dependence of the diffusion coefficient disappears when A = 2, because for a twoletter alphabet, every mutation implies a change in r.
For A > 2, there are mutations which do not change the Hamming distance and hence do not affect the diffusion process. This effect is reflected in the reduction of D in Eq. (10). Our continuous radial evolution equation (9) is somewhat reminiscent of the evolution equation in fitness space introduced by Tsimring, Levine, and Kessler (Tsimring et al., 1996) in a general population genetics context. However, with our concrete model for protein-DNA binding, we can work directly in genotype space, which will enable us to make explicit predictions on the behavior of the binding sites.

III. SELECTION THRESHOLD AND FUZZY MOTIFS
In this section, we use the evolutionary model (7) described in 'Model and Equations' to address the following questions: How large a selection pressure is needed for the maintenance of binding motifs? And can the fuzziness of the motifs be accounted for by the balance between mutation and selection? We will first provide an analytical solution to a simplified continuum model, and then show by numerical simulation that the qualitative features of the solution hold even for small system such as L = 10. We will compare these results to available data and discuss experimental ramifications.

Analytical results
Various properties of the continuum model specified by Eqs. (9-11) can be obtained exactly 3 . Here, we will present the results and discuss various qualitative features of the solution, in particular, the existence of a critical selection pressure for the maintenance of the binding motifs. Even though the continuum model is meaningful only for L ≫ r 0 ≫ 1, we will see from numerical simulation that the qualitative features are valid even for the more realistic parameter range where r 0 is not much larger than one, and L ∼ 10.
For the analytical study, we neglect the r-dependence of the diffusion coefficient (10) and the drift velocity (11), and use D = ν/2, v = ν. This is justified as long as r 0 ≪ L since as we will see, most of the interesting "action" of this system occurs around r = r 0 . Eq. (9) then reduces to the asymmetric "quantum well" problem well studied in the context of various statistical mechanics problems (Hatano and Nelson, 1997). [It differs from the DNA unzipping problem studied by Lubensky and Nelson (2000) only by an (unimportant) boundary condition at r = 0.] An explicit solution can be obtained by further approximating the Fermi function (4) by the Heavyside step function θ(r), such that the fitness landscape becomes This idealized form of the fitness function is known as truncation selection (Kondrashov, 1988). The solution of this simplified continuum problem is of the form n(r, t) = n 0 (r)e γt , where n 0 (r) is the stationary distribution associated with the largest growth rate γ. It is controlled by one dimensionless parameter, the effective selection pressure We have γ = 1 if α is below the critical value where η c is of order one and depends only weakly 4 on r 0 . In this regime, n 0 (r) is given by the continuum version of the (skewed) binomial distribution as if the fitness plateau at r < r 0 is not present. For α > α c , the solution is given in terms of the eigenvalue problem with the boundary condition y(0) = y ′ (0), where y(r) is the eigenfunction corresponding to the largest eigenvalue, λ( α), which must exceed 1.
[Thus, the precise definition of α c is λ( α c ) = 1.] In this regime, the growth rate becomes γ = 1 + (λ − 1)ν/2 ≥ φ 0 , and the stationary distribution is n 0 (r) = y(r)e r . The form of the latter can be straightforwardly analyzed for large α's (such that λ > 1). It is strongly peaked at r * < ∼ r 0 , indicating that the motifs are marginally conserved or maximally fuzzy 5 above the selection threshold.
A phase transition occurs at α = α c where the stationary distribution switches from being mostly confined in the region r < r 0 (localized phase) to the binomial distribution (delocalized phase) when λ approaches 1. This phase transition belongs to the same class of transitions as the one described by Eigen in the context of quasi-species evolution (Eigen, 1971;Eigen et al., 1989;Higgs, 1994). The critical selection pressure α c ∼ ν 0 ·L is recognized as the well-known form of the "error threshold". Note also the dependence of α c on r 0 as given in Eq. (14): α c decreases upon increasing r 0 , and since λ is a monotonously increasing function of α − α c , the effective growth rate γ will also increase. This implies that a wider fitness landscape has a smaller mutational load and a larger effective fitness, which is a known result, see, e.g. Schuster and Swetina (1988).
The "order parameter" of the phase transition is the average number of mismatches in the stationary state, r = r r n 0 (r)/ r n 0 (r). In the localized phase, r ≃ r 0 , while r = (A−1)L/A → ∞ in the delocalized phase. When α approaches α c from above, r diverges as 4 For 1 ≪ r 0 ≪ L, ηc is given to a good approximation by the the solution of the equation ηc = −r 0 · tan(ηc) and hence ηc ∈ [ π 2 , π] . 5 In the context of protein folding, it has been pointed out by R. Goldstein that the balance of mutation and selection may lead to maximal fuzziness in the space of amino acid sequences (Goldstein, 2001). Our results are similar in spirit, but more explicit due to the simplicity of the protein-DNA binding. indicating that this is a second-order phase transition 6 .

Numerical results
In order to test whether the behavior derived above for the simplified continuum model holds approximately also for the discrete model (7) with realistic parameters, we performed a number of numerical studies. We determined the steady-state distribution N 0 (r) of Eq. (7) over a range of values of α for two sets of parameters: (a) a nearly-continuum model, with L = 1000, r 0 = 30, and the step function landscape (12); and (b) the discrete model with L = 10, r 0 = 3, and the Fermi function landscape (8) with ε = 2k B T . Fig. 2 shows the stationary distribution N 0 (r) for the discrete model in the delocalized regime ( α = 0.5), in the localized regime ( α = 3.0), and in the crossover 7 region in between ( α = 1.4). We see that N 0 (r) is peaked slightly below r 0 = 3 in the localized regime, and becomes indistinguishable from the binomial distribution (15) in the delocalized phase. Note that the distribution is broad in the crossover regime, which is consistent with the finding of a continuous second-order 6 It should be noted that both the critical value αc and the divergence of r near αc are modified if one explicitly includes the time dependence of the fitness landscape. In particular, if we take the fitness to be φ(r, t) ∝ f (t)P (r) (f (t) = 1 in the ONstate and f (t) = 0 in the OFF-state), with a stochastic f (t), then the evolution dynamics becomes equivalent to the class of timedependent depinning problems studied by Lubensky and Nelson (2000), with the critical behavior r ∝ ( α − αc) −2 . 7 With a small system such as L = 10, the critical point of the phase transition becomes a crossover region, in which the behavior of the stationary state changes smoothly. To make the comparison more quantitative, we next examine the order parameter r . Fig. 3(a) shows r 0 / r plotted as a function of α, for the discrete model (triangles) and the nearly-continuum model (diamonds). The nearly-continuum model displays a sharp transition at 8 α c ≈ 1.6. The sharp transition becomes a pronounced crossover for the discrete model, but still with a relatively well-defined threshold α c . The r 0 -dependence of α c is 8 The critical value for the nearly-continuum model deviates somewhat from the analytical result (14). This deviation is due to the r-dependence of v and D, which we neglected in the analytical solution of the continuum model. plotted in Fig. 3(b) over the relevant interval 1 < r 0 < 5 (here, we have defined the threshold α c as the value of α where the derivative of r 0 / r ( α) is maximal). We see that it is relatively insensitive to the precise value of r 0 , with α c ∼ 1 as given qualitatively by the formula (14).

Viral and bacterial evolution
We expect the selection threshold described above to be detectable in evolution experiments with RNA viruses. The total mutation rate ν for the binding site for RNA viruses is of the order 10 −3 ∼ 10 −4 for a binding sequence of length L = 10. Assuming that the fitness gain of the virus in the ON-state (i.e., the factor δφ A in Eq. (2)) is of the order of 1% ∼ 10%, then the effective selection pressure α = 2f δφ A /ν on the viral regulatory sequence becomes of the order α c ∼ O(1) if the fractional exposure f to the ON-state is set at a few percent level. By varying f over the range of several percents, we expect that the phase transition should be observable. Moreover, the anomalous dependence (see footnote 4) of the selection threshold on the temporal variation f (t) should also be observable by applying controlled temporal changes to the environment. The stationary distribution N 0 (r) itself can be monitored in principle by sequencing a reasonable number (say 100) viral regulatory sequences after stationarity is reached.
A very different situation is expected for the evolution of bacteria or even DNA viruses. The total mutation rate ν is of the order 10 −8 for bacteria and 10 −6 for DNA viruses. Consequently, α is expected to be four orders of magnitude larger than α c for bacteria and two orders larger for DNA viruses. What is the behavior of the discrete model at such large values of α ? In Fig. 4, we show the position of the peak r * of the distribution N 0 (r) obtained numerically as a function of α on a logarithmic scale. For values of α exceeding ∼ 140, we find the peak is pushed down to r * = 0, contrary to the fuzziness depicted in Fig. 3. This behavior is obviously an artifact of the specific feature of the Fermi function landscape used in (8): for very large α's, there is an incentive for the distribution to move towards small r's due to the very slight increase in the value of P (r) for smaller r's. But it is unreasonable to expect that the simple relation between the binding probability P (r) and the fitness function φ(r) assumed in this study to hold down to very small differences in P (r). Aside from various kinetic effects of binding and temporal variations of the environment that we have neglected, fluctuations due to finite population size (e.g., genetic drift) simply do not allow for the population to resolve the very small differences in fitness due to the small differences in P (r); see the theory of nearlyneutral evolution (Ohta, 1992). Thus, φ(r) should be effectively r-independent for small r's. This can be implemented by replacing φ(r) by a constant value φ(r ′ 0 ) when φ(r ′ 0 − 1) − φ(r ′ 0 ) is below some resolution limit (set by the effective population size of the organism.) 9 For low mutation rate (or large α's), this amounts to replacing the fitness function by an infinte square well: The stationary distribution obtained in this case depends only on the width of the well r 0 , and is shown in Fig. 5 for r 0 = 3. Note that it is highly peaked at r 0 as expected. Hence the binding sequence is fuzzy even as α → ∞. However, it is different from simply truncating the binomial distribution (15) for r > r 0 due to the mutational load, i.e., a fraction of the population with r = r 0 will 9 This modification of the fitness function should in principle also be applied to the case of RNA virus evolution in the vicinity of the phase transition. However, it wouldn't make much of a difference there because the distribution would be already peaked away from small r. receive an additional deleterious mutation and not make it to the next generation.

Comparison to known sites
It is useful to compare the above solution to biological data. Unfortunately, polymorphisms in the same binding sequence across different strains of a bacteria species are not yet readily available. Instead, we assume that the different binding sequences (of the same transcription factor) located in different regulatory regions across the genome may be viewed as a sample of the stationary binding sequence distribution. This is clearly not the case if the selection pressure is small, since close to the selection threshold, even small differences in selection pressure experienced by the different binding sequences will produce different binding sequence distributions; see Fig. 2. But this should not be a concern for bacteria since α ≫ α c there. An obvious candidate is the binding sequences for the the well-known E. coli global regulator CRP (also known as the catabolite activator protein, CAP), which is activated under low cellular glucose level (Saier et al., 1996). There are over one hundred CRP sites in the E. coli genome. We take from the Regu-lonDB database (Salgado et al., 2000) a list of 28 sites which are biologically confirmed binding sites and appear only once in a given regulatory region. (The case of multiple binding sites is discussed in the following section.) The drawback of using CRP sites is that CRP is hardly ever the only regulator in a target regulatory region, and interaction with other transcription factors can complicate the situation. Fig. 6 shows the histogram of the number of mismatches of these CRP sequences from the consensus sequence TGTGA......TCACA. While it peaks at r * = 2 ∼ 3 similar to the corresponding distribution of the infinitemesa model in Fig. 5, it is clear that the distribution of the CRP sites is broader. The few outliers at r = 0, 1, 4, 5 may well be due to direct or subtle interaction with other factors which we have not considered in this simple model. The existence of nearly equal peaks at r = 2 and r = 3 is more perplexing: According to our model, the distribution should be peaked at the largest possible r. [For the L = 10 sequence, entropy favors r = 3 over r = 2 by a factor of eight.] One possible cause of the discrepancy may be the deviation of the real binding energy matrix ε i (S) from the two state model. For instance, suppose the chemical potential µ in the ON-state is such that r 0 = µ/ε is slightly below 3. Then according to the pure two-state model, the maximum number of allowed mismatches is 2. However, small deviations in the binding energies from ε will allow a maximum of 3 mismatches in a subset of the L positions, thereby producing a distribution peaked at both r = 2 and r = 3 as shown in Fig. 6. The actual stationary distribution of r can be easily computed numerically if the energy matrix is known. However at present, the authors know of no example of a transcription factor whose binding energy matrix is measured and a large number of binding sequences are available.
A very different explanation of the data in Fig. 6 is the differential selection of each of the CRP motifs as alluded already in "Introduction". Specifically, one can envision a situation where the single "ON-state" assumption we adopted is not valid, and instead the cell coordinates a graded response to cellular glucose shortage, requiring different operons to turn ON at different (activated) CRP concentrations. [In this case, our assumption that the fitness function (1) has a simple linear dependence on the binding probability obviously breaks down.] Within this scenario, the distribution in Fig. 6 is solely a result of the functional need of the cell, and its resemblance to the statistical distribution of Fig. 5 would be fortuituous. Distinguishing between the different plausibilities is important and can be done by either sequencing the CRP binding sites in a variety of related strains to accumulate statistics on polymorphism for each site, or perform site-directed mutagensis to specific binding motifs and directly measure the fitness function. In general, one may expect to find that the form of the fitness function depends on the biological function of the binding site. In particular, the form (1) seems more likely to be appropriate for the case of transcriptional repressors than for activators, since repressors need to have a binding probability close to one, in order to efficiently suppress transcription from the promoter, which is active in the absence of the repressor. In the case of activators, the promoter has a very low basal level of transcription and even an activator with a relatively low binding probability can lead to a large effect on the transcription level.

IV. MULTIPLE BINDING SITES
It is well known that regulatory binding motifs often occur in doublets or even higher multiplets. For instance, the regulatory regions of the E. coli genes crp, dadA, dsdXA, fixA, glpFK, glpTQ, lac, manX, nagE, and tsx are some of the many regions that contain two CRP binding sequences. Here, we extend our model to account for the possibility of multiple sites that bind the same protein and regulate the same promoter. We will pursue the question of whether we can interpret regulatory regions with multiplets as being under higher selective pressure for factor binding than regulatory regions with single binding sites.
Some factors (e.g., λ-repressor) bind cooperatively to binding sites, thereby effectively enhancing their DNA binding specificity. Cooperative factor binding can play an important and interesting role in transcription regulation (see, e.g., Ptashne (1992)), however it does so only for a fraction of the known multiplets, since many factors (such as CRP) have no binding domain for an attractive interaction between themselves. In the present study we exclude factor-factor interactions and explore possible selective advantages of multiple independent binding sites. This approach is similar in spirit to studies of gene duplication, which consider the evolution of multiple copies of the same gene (see, e.g., Wagner (2000)). One scenario could be that several bound transcription factors can simultaneously interact with polymerase to synergistically recruit (or repel) it more efficiently than a single factor would. For the case of CRP, this effect has been observed, and studied in detail, experimentally (Busby and Ebright, 1999;Langdon and Hochschild, 1999). An individual organism with a multiplet of binding sites for a factor then has a fitness advantage over one with a single binding site, if a strong activation (or repression) is beneficial for the biological function. Consequently, selection would favor multiplets over singlets. On the other hand, random mutations tend to destroy the binding motifs, so that an additional motif is maintained only when its contribution to the fitness is sufficiently high. In the following, we explore this scenario within our model.
Let us assume there are two binding sites in a certain regulatory region and ask whether they will be maintained by evolution. We begin by constructing a 'twosite fitness function' that makes the selection mechanism outlined above explicit. As in the previous sections on the single-site problem, we assume that the state of the bacterium/virus alternates between an ON-state, where factor binding leads to a fitness gain, and an OFF-state, where factor binding has a negative effect. In the ONstate, let the fitness gain due to factor binding be δφ A1 , δφ A2 , or δφ A12 , if a factor is bound to site 1 only, site 2 only, or both sites, respectively. Using the same arguments as for the single site case, the time-averaged fitness becomes where P 1 , P 2 denote the probabilities that a factor is bound to site 1, 2, which depend on the respective sequences (we neglect the possibility of cooperative binding at this point). The selection pressure, α, has again the form (2), with α = f δφ A1 , while the 'synergism factor' ω describes the fractional fitness advantage of two bound factors over just one, i.e. ω = δφ A12 /δφ A1 . In the remaining term, the dimensionless coefficient σ constitutes an 'asymmetry factor' equal to the relative fitness gain δφ A2 /δφ A1 (i.e. when σ = 1 one may distinguish between a more "important" and a less important site). Note that not only the selection pressure, but also ω and σ may vary between different regulatory regions, even when they are controlled by the same factor, since both depend on the location of the binding sites with respect to the promoter and on the sequence of the promoter itself (see the discussion at the end of this section).
As in the single site problem, we work in the two-state model approximation (see section 'Binding probability'), so that the binding probabilities P 1 and P 2 only depend on the number of mismatches r 1 and r 2 in the respective site, and take the form (4). When the selection pressure α is much larger than the mutation rate ν (as we typically expect in the case of bacterial evolution), we again invoke the argument that very small differences in the fitness function are hardly resolvable by finite populations, and therefore the fitness function should become neutral, i.e. r-independent, at small r 1 and r 2 . This effectively amounts to using step functions for the binding probabilities, i.e. P (r 1,2 ) = 1 for r 1,2 ≤ r 0 and P (r 1,2 ) = 0 for r 1,2 > r 0 . The two-site fitness function in (r 1 , r 2 )-space is then In order to simplify the discussion in the following, we will use the fitness function (20) over the whole range of α, since it yields a correct description at large α, and at small α, it produces no qualitative changes in the behavior of the stationary distribution compared to the smooth fitness function with P (r) of the form (4). It is straightforward to derive a two-site evolution equation analogous to Eq. (7), which describes the approach of the average distribution of mismatches N (r 1 , r 2 , t) (neglecting fluctuation effects due to finite population size) to its stationary state N 0 (r 1 , r 2 ). One obtains ∂ ∂t N (r 1 , r 2 , t) = φ(r 1 , r 2 )N (r 1 , r 2 , t) + + ν 0 ∆ r1 (r 1 +1)N (r 1 +1, r 2 , t) + − (A − 1)(L − r 1 )N (r 1 , r 2 , t) + ν 0 ∆ r2 (r 2 +1)N (r 1 , r 2 + 1, t) + − (A − 1)(L − r 2 )N (r 1 , r 2 , t) . (21) In the continuum limit, Eq. (21) becomes a twodimensional generalization of the (biased) diffusion equa- , v(r2)) in the continuum equation (22), while the shading corresponds to the fitness function (dark means high fitness; here we used ω = 2 for the purpose of illustration, and σ = 1) tion (9), ∂ ∂t n(r 1 , r 2 , t) = ϕ(r 1 , r 2 ) n(r 1 , r 2 , t) + + ∂ ∂r 1 D(r 1 ) ∂ ∂r 1 n(r 1 , r 2 , t) − v(r 1 )n(r 1 , r 2 , t) + + ∂ ∂r 2 D(r 2 ) ∂ ∂r 2 n(r 1 , r 2 , t) − v(r 2 )n(r 1 , r 2 , t) , (22) where D(r) and v(r) are still given by Eqs. (10) and (11) and n(r 1 , r 2 , t), ϕ(r 1 , r 2 ) are the continuum generalizations of N (r 1 , r 2 , t) and φ(r 1 , r 2 ), respectively. Fig. 7 illustrates the two-dimensional (biased) diffusion dynamics that emerges from Eq. (22). The fitness function has a high plateau or 'mesa' at r 1 , r 2 < r 0 , and two strips of lower fitness along the r 1 and r 2 axis. Hence selection tries to keep r 1 , r 2 < r 0 . Mutation, on the other hand, drives the distribution towards the average number of mismatches in a random binding sequence, r 1 = r 2 = (A − 1)L/A, as indicated by the arrows in Fig. 7. We are interested in the stationary distribution N 0 (r 1 , r 2 ) that arises as a balance between selection and mutation. Below we characterize the dependence of N 0 (r 1 , r 2 ) on the effective selection pressure α = 2α/ν and the synergism factor ω numerically by iterating Eq. (21). However, we first anticipate the qualitative behavior of N 0 (r 1 , r 2 ) using the understanding of the single site problem gained in the last section. Let us neglect a possible asymmetry between the sites for the moment, i.e. we set σ = 1. It is clear that if α is below a certain threshold value, no motif will be maintained, i.e. the peak of the stationary distribution will be close to r 1 = r 2 = (A − 1)L/A. On the other extreme, when α is very large, the distribution will certainly be localized on the high fitness mesa, corresponding to two conserved binding motifs. By analogy with the single site case, we would expect the distribution to be maximally fuzzy in this regime, and hence the peak of the stationary distribution to be close to r 1 = r 2 = r 0 . What happens when α takes on intermediate values? Upon increasing α, the peak of the stationary distribution may either pass directly from r 1 = r 2 = (A − 1)L/A to r 1 = r 2 = r 0 or go through a state with only one conserved motif (see Fig. 8). Intuitively, which of these "pathways" is taken, should depend on the value of ω: when ω is small, the selective advantage of two conserved motifs over one is small and therefore a much higher selection pressure is needed to stabilize two motifs against mutations than just one, i.e. upon increasing the selection pressure the system passes from 0 to 1 to 2 motifs. Conversely, when ω is very large, two motifs are actually stabilized at lower selection pressures than a single motif would be, and hence the system passes directly from 0 to 2 motifs. We can estimate the value ω c at which the system switches between the two pathways: when (ω − 1) ≪ 1, the 1motif phase exists in an intermediate range of α's, i.e. α c1 < α < α c2 , where the lower critical value for the transition from 0 to 1 motif is approximately the same as in the single-site problem, i.e. α c1 ∼ 1, and the upper critical value is α c2 ∼ (ω − 1) −1 , since the transition from 1 to 2 motifs may be regarded as another singlesite problem with α replaced by (ω − 1)α. The system switches between the two pathways when α c1 = α c2 , and hence ω c ≈ 2.
In the 1-motif phase, the selection of either motif, at site 1 or site 2, is equiprobable as long as σ = 1. Corre- spondingly, N 0 (r 1 , r 2 ) has two equal peaks as indicated in Fig. 8. When σ < 1 the peak associated with site 2 will disappear, but otherwise the qualitative behavior of the model remains the same. For simplicity, we will keep σ = 1 from here on.
In order to examine the qualitative picture outlined above and to render it more quantitative, we now characterize the behavior of N 0 (r 1 , r 2 ) numerically using the parameters tailored to CRP, i.e. L = 10 and r 0 = 3. To determine N 0 (r 1 , r 2 ), we evolve an arbitrary initial distribution N (r 1 , r 2 , t = 0) using Eq. (21) until the stationary state is reached. Fig. 9 displays three such stationary distributions, one each in the 0-motif, 1-motif, and 2-motif phase (here, we used α = 0.2, 5, and 50, together with ω = 1.1). Besides justifying the schematic sketch in Fig. 8, it shows that the distributions both in the 1-and 2-motif phase are peaked at the "edge" r 0 = 3 and are therefore maximally fuzzy as in the single site problem.
Next, we focus on the transitions between the three phases. In Fig. 10, the average total number of matches, i.e. 2L− r 1 − r 2 (here . . . denotes an average over the stationary distribution), is plotted against α, again with ω = 1.1 (solid line). [Note that in this figure the y-axis is reversed compared to Fig. 4.] We observe that the total number of matches rises quite sharply around α = 1 and α = 10. These positions agree with our estimates α c1 ∼ 1 and α c2 ∼ (ω − 1) −1 based on the qualitative discussion above. (Note that since we work with a small 'system size' of L = 10, the transitions, which are sharp in the limit L → ∞, appear only as smooth crossovers.) In order to show that the rises are indeed caused by the transitions from the 0-motif to the 1-motif phase, and from the 1-motif to the 2-motif phase, respectively, we also plotted the number of matches in each site at the peak (at one of the two peaks in the 1-motif phase) of the stationary distribution in Fig. 10 (circles and triangles). This illustrates the typical behavior of the individual sites: the first site jumps from 2 ∼ 3 matches to 7 matches at α ≈ 1 and the second site does the same at α ≈ 10. Hence, the motifs clearly appear in a step-like fashion as the selection pressure is increased.
In evolution experiments with RNA viruses, this twofold transition should be directly observable (if a suitable operon can be identified where the fitness function (19) is applicable), since according to our estimates above, α for these systems can be tuned over a range of 1 ∼ 100 by varying the fractional exposure to different environmental conditions. On the basis of our model, one would expect, for instance, that one of the sites in a doublet disappears in the course of an evolution experiment, if the selection pressure is sufficiently lowered by reducing the exposure to the environment where binding is beneficial. When the exposure is reduced to zero, both regulatory sites and the gene coding for the transcription factor (if not required for other mechanisms) will be lost.
To complete our characterization of the model behavior, we map out the entire phase diagram in the ( α, ω) parameter-space. The result is shown in Fig. 11, where α c1 and α c2 are plotted as a function of ω. Since L = 10 in the present case and the phase boundary is well-defined only in the limit L → ∞, the curves α c1 (ω) and α c2 (ω) are really only crossover lines whose precise location is slightly dependent on their definition (in the figure they are represented by dashed lines in order to indicate this fact). [Here, we defined α c1 and α c2 as the value of α where the peak of the stationary distribution first reaches 7 matches in the respective site.] We see that the phase boundaries join at ω ≈ 2 as expected from our estimate above. When ω → 1 from above, the upper boundary, α c2 , diverges as (ω − 1) −1 . This implies that at ω = 1, the two-motif phase is unreachable, regardless of the value of α. Finally, we note that Hermisson et al. (2001) observe a somewhat similar two-fold transition between three phases in a situation where two different mutation processes with distinct mutation rates are taken into account. Let us now return to the scenario outlined at the beginning of this section, and discuss whether we can interpret the regulatory regions with multiplets as being under higher selective pressure for factor binding than the ones with single binding sites. Our study of the two-site problem would justify this conclusion, if (a) the values of ω and σ were very similar for all regulatory regions, and (b) the effective selection pressure were typically on the same order of magnitude as α c2 . Then, we could tentatively associate singlets with an α below α c2 , and multiplets with an α above that threshold. However, both conditions are not likely to be fulfilled generically in bacteria. First, the values of ω and σ not only depend on the sequence of the promoter and the distance of the binding sites from the promoter, but also on the level of gene expression that is beneficial for the cell. For example, genes that code for proteins which are not needed in large amounts, even in the environmental condition where their expression is favored, do not require a large activation, and hence the effect of a second binding site could even be detrimental, i.e. ω < 1. Second, α should typically be on the order of 10 4 (see above), and hence, as long as ω is only slightly larger than one (ω > 1.0001), our model would always predict multiple binding sites for bacterial transcription factors. Therefore, within our model, whether one or two motifs are maintained depends almost exclusively on the value of ω, i.e. ω ≤ 1 leads to one motif and ω > 1 to two motifs. However, there may be cases where maintaining subtle differences in the temporal ordering of turning on/off different operons would give rise to a very small fitness advantage, e.g. flagella assembly and SOS-respnse in E. coli (see recent results by U. Alon, submitted for publication). In such cases, the system may respond by keeping one or two motifs according to the theory we presented. And of course there is also the situation of RNA viruses described above where the twofold transition as depicted in Fig. 10 could in principle be directly observed.

V. DISCUSSION AND CONCLUSIONS
The fuzziness of regulatory binding motifs is a widely observed phenomenon. The present investigation has shown that the entropic advantage of introducing mismatches from the best binding sequence is sufficient to produce motifs that are maximally fuzzy while still retaining the capability of factor binding. Nevertheless, we cannot exclude that the fuzziness actually bears a selective advantage (in the language of population genetics, this would correspond to a stabilizing selective pressure). The alternative scenario given for the fuzziness of the CRP sites is an explicit example of the latter. It would be very interesting to address this question experimentally by directly measuring the fitness of a bacterium or virus as a function of the sequence of its binding sites: starting with a wild-type binding site that has several mismatches, what is the effect on the fitness, when the number of mismatches is reduced by site-directed mutagenesis? Does the fitness remain unchanged or is it reduced? Besides answering the question raised above, experiments of that type could also lead to important conclusions on the evolution of genetic regulation.
Another important result of our study is the phase transition associated with the maintenance of motifs. Our general conclusion is that the selection pressure on a single binding motif needs to surpass a threshold value of approximately ν 0 L/2 to guarantee maintenance, while the threshold for a second site (for the same factor and in the same regulatory region) is larger by a factor (ω−1) −1 , where ω is given by the ratio of the fitness of the organism with two sites over the fitness of the organism with one site. As pointed out above, this prediction could be tested experimentally by evolving RNA viruses in a fluctuating environment and varying the fractional exposure to the environment where factor binding is beneficial. In this case there would be no need to do site-directed mutagenesis, since the transition could be observed directly by sequencing.
Our model makes quantitative predictions on the statistics of polymorphisms in binding sites. To test these, it will not suffice to sequence a particular binding site in many different isolates from a single, large (N ν 0 ≫ 1) laboratory population, since this population originates from a small, genetically homogeneous population and it takes a time on the order of ν −1 0 to equilibrate the distribution of mismatches in a binding site. Instead, sequencing the same binding sites in several different strains should yield the desired data. Besides allowing a comparison to our model, detailed information on polymorphisms in binding sequences would also make it possible to address a number of interesting questions, e.g. how does the selection pressure on binding sites compare with the selection pressure on coding sequences 10 ? Or, can one identify compensating mutations between promoter and binding site sequences, e.g. could a mutation that weakens the promoter be compensated by a mutation that increases the affinity of an activator to its operator site?
We conclude that the evolution of transcription factor binding sites is a problem that is well suited to establishing a link between the detailed molecular biophysics of a system and its evolutionary dynamics. The theory presented in the present article is meant as a first step, with the hope of stimulating future experiments in this direction. There are many directions for the improvement of the model and the analysis. One important question is the validity of the mean-field analysis described here. Is the finite population size effect important here and how would it change the motif statistics within our model? One can also investigate more elaborate models including, for instance, the effect of time-dependent environment, the coupled evolution of the polymerase and transcription factor the binding sites, and the interaction among different factors.