Adaptation, fitness landscape learning and fast evolution

We consider evolution of a large population, where fitness of each organism is defined by many phenotypical traits. These traits result from expression of many genes. Under some assumptions on fitness we prove that such model organisms are capable, to some extent, to recognize the fitness landscape. That fitness landscape learning sharply reduces the number of mutations needed for adaptation. Moreover, this learning increases phenotype robustness with respect to mutations, i.e., canalizes the phenotype. We show that learning and canalization work only when evolution is gradual. Organisms can be adapted to many constraints associated with a hard environment, if that environment becomes harder step by step. Our results explain why evolution can involve genetic changes of a relatively large effect and why the total number of changes are surprisingly small.


Introduction
A central idea of modern biology is that evolution proceeds by mutation and selection. This process may be represented as a walk in a fitness landscape leading to fitness increase and slow adaptation 1 . According to classical ideas this walk can be considered a sequence of small random steps with small phenotypic effects. Nevertheless, there is a limited amount of experimental support for this idea 2 and some experimental evidence that evolution can involve genetic changes of a relatively large effect and that the total number of changes are surprisingly small 3 . Another intriguing fact is that organisms are capable of making adaptive predictions of environmental changes 4 .
To explain those facts new evolutionary concepts have been suggested (see the review by 5 and references therein). The main idea is that a population can "learn" (recognize) fitness landscapes 5-7 . The approach developed in these works is a generalization of ideas from machine learning in which learning (regression to data) is viewed as selection and generalization (interpolation or extrapolation) is viewed as adaptation.
A mathematical basis for investigation of evolution learning problems has been developed by 8. However, this work uses a simplified model, where organisms are represented as Boolean circuits seeking an "ideal answer" to environmental challenges. These circuits involve N g Boolean variables that can be interpreted as genes, and the ideal circuit answer maximizes the fitness. A similar model was studied numerically by 7 to confirm the theory of "facilitated variation" explaining the appearance of genetic variations which can lead to large phenotypic ones. In the work by 9 a theory of the evolution of these Boolean circuits was advanced. It was shown that, under some conditionsweak selection, see 10-a polynomially large population over polynomially many generations (polynomial in N g ) will end up almost surely consisting exclusively of assignments, which satisfy all constraints. This theorem can shed light on the problem of the evolution of complex adaptations since that satisfiability problem can be considered as a rough mathematical model of adaptation to many constraints.
In 6 it is shown that, in the regime of weak selection, population evolution can be described by the multiplicative weight update algorithm (MWUA), which is a powerful tool, well known in theoretical computer science and a generalization of such famous algorithms as Adaboost and others 11 . Note that in 6 infinitely large populations are investigated whereas the results of 9 hold only for finite populations and take into account genetic drift.
Evolution Computation(EC) problems are considered recently by many papers 12-16 mainly for artificial test fitness functions like OneMax or LeadingOnes (for an overwiew of EC problems, see 17).
In this paper, we investigate adaptation and the fitness landscape learning problem for more realistic fitness function This fitness can model adaptation for insects and connected with a fundamental hard combinatorial problem: K-SAT.
The main results can be outlined as follows. We show that, in a fixed environment, genes can serve as learners in the machine learning sense. Indeed, if an organism has survived for a long period, this fact alone constitutes important information, which can be used. The biological interpretation of this fact is simple: if a population is large enough and mutations are sufficiently rare, deleterious mutations are eliminated by purifying selection. Hence, those non-neutral mutant alleles which have become fixed in natural populations will, with probability close to 1, be adaptive and cause a positive increment of fitness (see Theorem 3.1 and Theorem 3.2 in Subsection 3.1 and Subsection 3.2). We obtain mathematical results, which allows us to estimate the reduction of mutation number due to that learning landscape procedure. Learning can sharply reduce the number of mutations needed to form a phenotypic trait useful for adaptation that is consistent with experimental data mentioned above (see 3).
Another important result is as follows. We estimate the accuracy of fundamental Nagylaki equations 6,10 for a realistic population model, where the population size is bounded and a non-zero mutation rate is taken into account (in the case of asexual reproduction). Those accuracy estimates are fulfilled for all possible values of mutation rates and population sizes.

Model
In this section, we describe our model and mathematical approach.

Genome
We assume that the genotype can be described by Boolean strings of length N g , where N g is the number of genes. Then where s i = 1 means that gene i is activated (switched on) and s i = 0 means that it is repressed (switched off). Correspondence

Amendments from Version 1
To respond to the reviewers' comments, in this revision, we extended the discussion. For the key point that the network modifies the thresholds by a simply feedback mechanism we point to regulation via enhancers.
In order to explain in more detail the regulation we added a new Section 2.5 on "Gene regulation networks". Two new figures ( Figure 1 and Figure 2) show results of numerical simulations based on the ``strong selection weak mutation'' (SSWM) algorithm.
In the discussion of Theorem 3.2 we added a paragraph on rare mutants and giving some intuitions towards the proof of Theorem 3.2 as it is based on estimates of the accuracy of the Nagylaki equations. For the main Lemma 4.5 for the proof of Theorem 3.2. we add a transparent interpretation.
Additional references pointed out by the reviewers are incorporated into version 2 of the paper.

Any further responses from the reviewers can be found at the end of the article
between Boolean hypercube and genotypes is considered for example in 12.

Phenotypic traits
Although phenotype is controlled by genes, it is also influenced by environmental conditions and various epigenetic processes.
In this paper, we suppose that phenotypic traits are controlled by genotype only. We consider levels f j of expressions of those traits as real variables in the interval (0, 1). Then the vector f = (f 1 , . . . , f N b ) can be considered to represent the organismal phenotype. We suppose that where f j ∈ (0, 1) is a real valued function of the Boolean string s, the genotype.
Only a part of s i is involved in f j . Namely, for each j we have a set of indices K j = {i 1 , i 2 , . . . , i nj } such that f j depends on s i with i ∈ K j , so that where i l ∈ K j and n j is the number of genes involved in the control of the trait expression.
The representation of phenotype by the quantities f j is suggestive of quantitative traits because the f j are real valued. The limiting values of 0 or 1 suggest another interpretation, however, in terms of cell type. Multicellular organisms consist of cells of different types. One can suppose that the organismal phenotype is defined completely by the corresponding cell pattern. The cell type j is determined by morphogenes, which can be identified as gene products or signaling molecules that can change cell type (or genes that code for signaling molecules that can determine cell types or cell-cell interactions and then finally the cell pattern). The morphogene activity is defined by (2.2).
We further suppose: Assumption M. Assume activities f j have the following properties.
The sets K j are independent uniformly random subsets of We denote the total number of genes involved in regulation of all f j by N r , where Assumption M implies that the genetic control of the phenotype is organized, in a sense, randomly, and that only a portion of the full set of genes controls phenotypic traits. That modularity of gene control is well known from experimental data (see 18,19) and for evolution computation problems it was studied, for example, in 16.
Consider an example, where the assumption M holds, where we have a saturated expression, inspired by earlier work 20,21 . Let and w ij , h j are some coefficients (their meaning will be explained below). As an example, we can take σ (S) = (1 + exp(−bS)) −1 , where b > 0 is a sharpness parameter. Note that for large b this sigmoidal function tends to the step function and for b = +∞ our model becomes a Boolean one. The parameters h j defines thresholds for trait expression 20 . The relation (2.4) can be interpreted as a simple mathematical model for quantitative trait locus (QTL) action.
To understand the role of h j consider a trait f j and suppose that for a well adapted organism f j ≈ 1. Let, for simplicity, w ij take the values 1, 0, or −1. Then the parameter h j defines how many genes involved in the control of the f j expression should be activators and how many should be repressors. Let the numbers of activator and repressor genes be j n ± , respectively.
One can suppose that h j describes a direct influence of environment on phenotype, such as stress, that can exert epigenetic effects. In Section 2.6 using data from 22 we will show that the model defined by (2.4) are capable to describe main topological characteristics of really observed fitness functions in the case of mimicry, camouflage and thermoregulation for insects.
Let us introduce the matrix W of size N b × N g with the entries w ij . The coefficients w ji determine the effects of terminal differentiation genes (see 23), and hence encodes the genotype-phenotype map. We assume that the coefficients w ji are random, with the probability that w ji > 0 or that w ji < 0 is β/2N, where β > 0 is a parameter. This quantity β << N defines a genetic redundancy, i.e., averaged numbers of genes that control a trait. Note that then large β  1 one has n j < 4β with the probability Pr β , which is exponentially close to 1: Pr β > 1 -exp(−0.1β), thus, the number n j are bounded.

Fitness
We know little about the details of how fitness relates to the phenotype of multicellular organisms, and for that reason classic neo-Darwinian theory takes fitness to be a function of genotype. Some models which take account of epistasis have been proposed 24 . The random field models assign fitness values to genotypes independently from a fixed probability distribution. They are close to mutation selection models introduced by 25, and can be named House of Cards (HoC) model. The best known model of this kind is the NK model introduced by Kauffman and Weinberger 26 , where each locus interacts with K other loci. Rough Mount Fuji (RMF) models are obtained by combining a random HoC landscape with an additive landscape models 27 . In evolution computations (EC) some artificial fitness models were used, for example OneMax and Leading Ones to test evolution algorithms, see for example 15 .
In this work, we use the classical approach of R. Fisher by introducing an explicit representation of phenotype, f, and allow it to determine fitness through interaction with an environment b. That is, we assume that the phenotype is completely determined by the phenotype trait expression, and thus the fitness depends on the genotype s via f j .
We express the relative fitness F and its dependence on environment b via an auxiliary function W via the relation ( , ) exp ( ( , )), where K F is a positive constant and b = (b 1 ,..., b N b ) is a vector consisting of coefficients b j , respectively. Below we refer to W as a fitness potential, and we assume that Sometimes, if the parameter b is fixed, we shall omit the corresponding argument in notation for W and F.
We consider fitness as a numerical measure of interactions between the phenotype and an environment. For a fixed environment, this idea gives us the fitness of classical population genetics. A part of the fitness, however, depends on the organism developing properly and for now we represent it as independent of the environment, although we are aware that this is not always the case. Note that some coefficients b j may be negative and others may be positive, and that the model ( where f * is a target pattern corresponding to an insect to mimic. Colour patterns can be also described by classical RGB formalism. The representation of the fitness as a sum of terms (2.7) is of course a rough approximation; however if assumption M holds that representation is consistent with important observed facts. First, mutations have been identified that alter one part of the pigment pattern without affecting any other. This independence of different pattern parts can be explained by the modular organization of the genetic regulation that controls pigmentation. In the course of evolution, different aspects of the pigment pattern have clearly evolved independently of each other 18 . Second, the topology of the fitness landscapes was studied in 22  In Section 2.6 we will show that the fitness model defined by (2.4) and (2.7) have those topological properties.

Population dynamics model
For simplicity, we consider populations with asexual reproduction. (Although a part of the results remain valid for sexual reproduction, as we discuss at the end of this subsection). We choose initial genotypes randomly from a gene pool and assign them to organisms. This choice is invariant with respect to the population member, i.e,. the probability to assign a given genotype s to a member of the population does not depend on that member.
In each generation, there are N pop (t) individuals, the genome of each of which is denoted by s(t), where t = 0, 1, 2,... stands for the evolution step number). Following the classical Wright-Fisher ideas, we suppose that generations do not overlap. In each generation (i.e., for each t), the following three steps are performed: 1. Each individual s at each evolution step can mutate with probability p mut per gene; 2. At evolution step t each individual with a genotype s produces k progeny, where k is a random non-negative integer, distributed according to the Poisson law where q = F(s) is the fitness of that individual; 3. To take into account ecological restrictions on the population size, we introduce the maximal population size N popmax . If N′(t) > N popmax , where N′(t) is the number of progeny produced by the population at step t, we kill randomly selected individuals in a population-dependent manner. The probability of the death of an individual is given by p kill (N′) = 1 − (N popmax /N′(t)). If N′(t) ≤ N popmax , we do nothing. We refer to this as the "massacre procedure." Conditions 1 and 2 imply that mutations in the genotypes create a new genetic pool and then a new round of selection starts. Condition 3 expresses the fundamental ecological limitation that all environments can only support populations of a limited size. If N popmax  1 then by (2.8) and the Central Limit Theorem one can show, under some additional conditions, (see Section 4) that fluctuations of the population size are small, and thus the population is ecologically stable and N pop (t) ≈ N popmax .
In the limit case of infinitely large populations we will write the discrete dynamical equation for the time evolution of the frequency X(s, t) of the genotype s in the population as 1 0 where F(t) is the average fitness of the population at the moment t defined by The equations (2.9) do not take mutations into account. They only describe changes in the genotype frequencies because of selection at the t-th time step. The same equations govern evolution in the case of sexual reproduction in the limit of weak selection 6,10 . Note that for an evolution defined by (2.9), the average fitness F(t) defined by (2.10) satisfies Fisher's theorem, so that this function increases at each time step t: F(t + 1) ≥ F(t).

Gene regulation network
In this section, we follow ideas of the classic paper 29 : the model should include a regulatory network, which evolves itself.
Regulatory genes as well as environmental factors, such as temperature, can influence the trait expression. This effect can be realized via thresholds h j (we shall describe it below), or via a regulation of coefficients w ij (see 30). In fact, these approaches are similar for sharp sigmoidal functions σ that are close to step functions, as can be shown by ideas from 31. Consider the expression involved in relation (2.4). Suppose following 31 that w ij take the values γ, 0 or -γ, where γ is a parameter, which can be regulated. Let h j  γ be fixed.
Assume that at an evolution step we have S i > d S , where d S > 0 is a parameter, which is more than γ. According to our Theorems this fact indicates, that for a well adapted organism the trait f j (s) ≈ 1 and the corresponding coefficient b j = 1. On the contrary, if S i < -d S , then one can expect that b j = -1 and f j (s) must be 0.
In the both cases, we can regulate the trait expression by a feedback so that whenever S i attains a critical level, i.e., a trait is well expressed, then γ should be increased. In our numerical simulations we use an alternative model, where we change h j . The alternative model, which is used in our numerical simulations, can be described as follows. We suppose that depending on activity of some regulatory genes or proteins (such as Hsp90), the threshold value can take three values ( ) where D h > 1 is a large parameter that defines the number of genes involved in trait specification (see Subsection 2.2). Thus, h j can take large negative or positive values, and also a neutral value close to 0. The feedback can be described as follows: In our simulations, each ∆T evolution steps we modify h i from 0 to -D h or D h for the trait with the maximal value i S .
Such a regulation produces a good adaptation even when the number of genes is essentially less than the number of the traits.
The plot in Figure 3 shows a difference between an evolution without any regulation and with the regulation by (2.12) described above.
However, the evolution of gene regulation via h i has an advantage: it makes phenotype more robust. In fact, the traits with large i h are non-sensitive with respect to mutations. The effect produced by this robustness is shown on Figure 1.
This regulation is even more effective, if we consider the modular evolution, following the recent paper 32 . Indeed, biological systems are characterized by a high degree of modularity. This modularity allows biological systems to vary only in a small subset of traits at each evolution round. Figure 2 shows an effect of this modularity. We consider a toy example, where a system with only 10 genes should be adapted to 200 constraints. Without regulation, we have no chances to make an adaptation (only about 10 traits are correctly adapted, see the green curve). It is a consequence of a formidable pleiotropy (20 traits on 1 gene). However, if at each evolution stage an organism should be adapted to 4 traits whereas the remaining ones are made robust by high i h , then already 66 traits are correctly expressed after 20000 evolution steps. Note that the learning plays a key role in the regulation. In fact, the sign of the regulation threshold h i depends on the sign of the corresponding coefficient b i , and the knowledge of that sign gives us important information on the fitness landscape.

Adaptation as a hard combinatorial problem
Adaptation (i.e., maximization of fitness in a changing environment) is a very hard problem since over evolutionary history we observe the coevolution of many traits accompanied by changes in many genes. In its general context, this is a problem in the theory of macroevolution, which in general requires the integration of population genetics and developmental biology for its full understanding. There are two key components of this problem. Figure 1. This graph illustrates a difference between the adaptation for evolution without evolution of gene regulation via threshold (the red curve) and with that evolution by 2.12 (the green curve). The parameters are N g = 50, M = 300, the mutation rate p mut = 0.01, γ = 1 and K = 4. Non-zero coefficients w ij are random numbers distributed according to the standard normal law. The initial genome is a random binary string, where each value is 0 or 1 with probability 1/2. The coefficients b i are either 1 or -1, where the probability of 1 is p b = 0.8.
First, development is itself a dynamical process operating over time. Second, there is a combinatorial component of development wherein different combinations of gene must be expressed in different cell types. This combinatorial aspect of the problem means that straightforward theoretical methods of considering the relationship between gene expression and a changing environment that have been very successful in single celled organisms 33 cannot be applied to metazoa. In this work, for the sake of tractability, we focus on the combinatorial aspect of the problem and neglect developmental dynamics. Even at the highly simplified level of our model, adaptation is a hard computational problem, as we now demonstrate.
Consider the case, where f j are defined by relations (2.4) and assume that i) σ is the step function; As a consequence of the second assumption, F attains its maximum for f 1 = 1, f 2 = 1,..., f N b = 1. Let us show that, even in this particular case, the problem of the fitness maximization with respect to s is very complex. In fact, for a choice of h j it reduces to the famous NP-complete problem, so-called K-SAT, which has received a great deal of attention in the last few decades (see 34-39). The K-SAT can be formulated as follows. The problem is to test whether one can satisfy all of the clauses by an assignment of Boolean variables.
Cook and Levin 35,34 have shown that the K-SAT problem is NPcomplete and therefore in general it is not feasible in a reasonable running time. In subsequent studies-for instance, by 36-it was shown that K-SAT of a random structure is feasible under the condition that The set C K of solutions of random K-SAT has a nontrivial structure depending on parameter α = N b /N g 37,39 . For sufficiently small α < α g (K), where α g (K) ≈ 2 K log(K)/K is some critical parameter, the set C K forms a giant cluster, where nearest solutions are connected by a single flip and one can go from a solution to another by a sequence of single flips (pointed mutations) 39 . For α ∈ (α g , α d ), where α d (k) < α c is another critical value, solutions form a set of disconnected clusters. The local search algorithms do not work in the domain α > α g . Figure 2. This graph shows that the modular evolution of gene regulation allows adaptation with a few genes, in that case N g = 10 and the number of traits M = 200. This means that we have a very big pleiotropy. Evolution proceeds in 20000 steps. The green curve corresponds to random walk with fixed small h without any evolution of gene regulation. We observe that with 10 genes 9 traits are correctly expressed. If evolution goes into 50 rounds then 66 traits are correctly adapted, and if evolution goes in 190 rounds then 129 traits are correctly adapted (the red curve). In that last case, at each step we make adaptation to at most one trait. Parameters are as follows. The mutation rate p mut = 0.01 and K = 5. Non-zero coefficients w ij are random numbers distributed according to the standard normal law. The initial genome is a random binary string, where each value is 0 or 1 with probability 1/2. The coefficients b i are either 1 or -1, where the probability of 1 is p b = 0.8.Thus such mutations produce non-viable organisms.
Probably, for evolution context K-SAT was applied first in 40, where it was used for an investigation of speciation problem.
To see the connection of our model with K-SAT, consider equation (2.4) supposing that w ij ∈ {1, 0, −1} and h j = −C j + 0.5, where C j is the number of negative w ji in the sum 1 . maximize the fitness, we must assign s j such that all disjunctions will be satisfied. If we fix the number n j of the literals participating in each disjunction (clause) and set n j = K, this assignment problem is precisely the K-SAT problem formulated above.
Reduction to the K-SAT problem is a transparent way of representing the idea that multiple constraints need to be satisfied. The number K defines the gene redundancy and the probability of gene pleiotropy. Remind that pleiotropy occurs when one gene influences two or more seemingly unrelated phenotypic traits. The threshold h j and K define the number of genes which need be flipped in order to attain a high expression of the trait f j .
Note that gene pleiotropy is a fundamental characteristics 41 , which is studied for real organisms only recently (see 19). We can compare experimental observations and consequence of model (2.4), which is a generalisation of K-SAT (compare plot Figure 3 and plots on Figure 1 in 19). So, we can fit our model parameters using real data. Moreover, we can check validity of our model by the following arguments.
We note that, in the case of giant cluster formation, the topological properties of the solution set C K , mentioned above, outline the properties of really observed fitness landscapes 22 : existence of many peaks, valleys and ridges connecting peaks. Namely, existence of many solutions of K-SAT, when a giant cluster exists, means that the landscape has a number of peaks separated by valleys. On the other hand, connectance of solutions within the giant cluster can be interpreted that there exist ridges that connect peaks.
Note that there are important differences between K-SAT in Theoretical Computer Science and fitness maximization problems. First, the signs of b j are unknown for real biological situations since the fitness landscape is unknown. Second, our adaptation problem involves the threshold parameters h j (see (2.4)). In contrast to K-SAT, in our case the Boolean circuit is plastic, because the h j are not fixed.
If the b j are unknown, the adaptation (fitness maximization) problem becomes even harder because we do not know the function to optimize. Therefore, many algorithms for K-SAT are useless for biological adaptation problems. Below we will nonetheless obtain some analytical results based on the assumption that b j are random.

Main theorems
The subsequent material is organized as follows. First we formulate a result on regulation mechanism power. Furthermore, we prove two fitness landscape learning theorems.

Fitness landscape learning theorems
For simplicity, we consider asexual reproduction. To obtain similar results for sexual reproduction, one can consider a weak selection regime and use the results of 10, where eq. (2.9) are derived.
Let us introduce two sets of indices I + and I − , such that I + ∪ I − = {1,..., N b }. We refer to these sets in the sequel as positive and negative sets, respectively. We have The biological interpretation of that definition is transparent: the expression of the traits f j with j ∈ I + increases the fitness and for j ∈ Iexpression of the trait decreases the fitness. We formulate two theorems on fitness landscape learning. First we consider the case of infinitely large populations. we have j ∈ I + . If f j (s) > f j (s ), then j ∈ I -.
Before proving this, let us make some comments. The biological meaning of the theorem is simple: for simple fitness models, where unknown parameters b j are involved in a linear way, in the limit of infinitely large populations fitness landscape learning is possible.
Moreover, note that we do not make any specific assumptions about the nature of mutation, but only that all genetic variation between s and s are contained in a single regulatory set K j .
The assertion of Theorem 3.1 is not valid if the set Diff(s, s ) belongs to a union of different regulation sets K j , j = j 1 , . . . , j p with p > 1. This effect of belonging to different sets K j is pleiotropy in gene regulation. Note that if N b  N g then the pleiotropy probability is small for large genome lengths N g . On the contrary, if N b  N g then assumption II is invalid.
Assumption II looks natural if when we deal with point mutations. In fact, if s is obtained from s by a single point mutation then condition (3.4) always holds for some j. For small mutation rates the probability of two point mutations is essentially below than the probability of a single mutation.
To conclude let us note that Theorem gives a rough estimate for the learning time T c : Proof. The main idea is simple. Negative mutations lead to elimination of mutant genotypes from the population, and the corresponding frequencies become, for large times, exponentially small.
Assume that (3.7) holds. Let j ∈ I − , and thus b j < 0. Consider the quantity According to assumption II Assumption III entails that Relations (2.6) and (3.10) imply By (2.9) and the last inequality we find that for T > T 1 (3.11) Consider inequality (3.11) for T = T 1 + T c . Let us note that in the relation Q(T 1 ) = X(s, T 1 )/X(s , T 1 ) the numerator is p 0 whereas the denominator ≤ 1. Thus, Q(T 1 ) ≥ p 0 . The same arguments show that Q(T 1 + T c ) ≤ 1/p 1 . Therefore, by (3.11) one obtains that This inequality leads to a contradiction for T c satisfying (3.6), thus completing the proof.

The case of finite populations
Theorem 3.1 can be extended to the case of finite populations and non-zero mutation rates. is small. To formulate this generalization, we need an additional assumption about the fitness function. Suppose that Condition (3.13) means that each individual gives birth to at least c F and at most C F descendants, where those bounds do not depend on the population size and evolution step. It is interesting to compare Theorem 3.1 and Theorem 3.2. The previous one asserts that for infinite populations the probability of the event j ∈ I − is zero whereas the second one claims that this probability becomes exponentially small as the population size increases.
This theorem also shows that evolution can make a statistical test checking the hypothesis Hthat j ∈ Iagainst the hypothesis H + that j ∈ I + . Suppose that His true. Let V be the event that the frequency X(s , t) of the genotype s in the population is larger than p 1 within a sufficiently large time T c . According to estimate (3.18), the probability of the event V is so small that it is almost unbelievable. Therefore, the hypothesis Hshould be rejected. We will refer T c as the checking time.
Rare mutants. In this Theorem we assume that the frequencies p 0 and p 1 of genotypes (wild and mutant) are fixed and our estimate is valid as p mut → 0. I.e., we do not consider mutants with a very small frequencies (fractions). Of course, a large population always contains a small number of such mutants. In numerical simulations we assume that evolution is successful and population is perfectly adapted, if, say, 95 or 99 percents of population members have the maximal fitness.
Ideas for the proof. The main idea is the same as that for the previous theorem: we compare the frequencies of the organisms with the genotype s and the organisms with the genotype s.
However, the proof includes a number of technical details connected with estimates of mutation effects and fluctuations. The formal proof can be found in Section 4. It is based on estimates of the accuracy of the Nagylaki equations (2.9). The main Lemma 4.5 for the proof of Theorem 3.2. admits a transparent interpretation. We show that fraction X(s, t) of genotype s evolves in time in such a way that the estimates mut popmax is a fitness of genotype s, F is average population fitness, and r(p mut , N pop )) are small corrections, which converge to zero uniformly in s as the mutation rate p mut → 0 and the population size N popmax → ∞. This means that in the limit p mut → 0, N popmax → ∞ we have equation (2.9). The main problem with the application of Theorem 3.2 is how it allows to perform fitness landscape learning. It can be done by a regulation, as is detailed in the following section.

Proof of theorems
Let us prove Theorem 3.2.

Main tools and auxiliary Lemmas
Let us introduce notation and make some preliminary remarks. Remind that we denote by N(s, t) the number of the population members with the genotype s at the moment t. Let X(t) be the set of all population members at the moment t. For each x ∈ X(t) let us denote by N′(x, t) the number of progeny born by the individual x at the moment t before the massacre (see point 3 of model from Subsection 2.4). Let s g (x) be the genotype of x. Then, according to (2.8), the mean of N′(x, t) is where EX denotes the expected value of X. By

( )
, N s t we denote the number of all progeny born by individuals with the genotype s at the moment t before the massacre. Since all progeny are produced independently and randomly, the previous relation gives Our main analytical tools are the Chernoff bounds and the Hoeffding inequalities. We also use the Markov inequality: for a positive random quantity X and a > 0 one has We minimize f with respect to λ and obtain (4.6). To derive (4.7), we use ∑ of the mutually independent random quantities. According to (4.2), the average EN' (x, t) is F(s g (x)). Therefore, Then one has (4.28) Let us define the event (4. 29) where the interval J ε (t) is defined by (4.16) and ε  is defined by (4.27). By (4.4) we have 2 2 Pr (Not ( )) Pr ( | ) Now we apply the Hoeffding inequality (4.15). For each ε 2 > 0 we obtain Therefore, Moreover, by Lemma 4.3 Inequalities (4.30), (4.32) and (4.33) prove (4.25).
The following lemma, in particular, allows us to obtain equations (2.9) and (2.10) in the limit of infinite populations and for small mutation probabilities.
Recall that ( , ). N s t denotes the number of non-mutated progeny generated by the individuals with the genotype s before the massacre. Let N sur (s, t) be the number of those progeny that survived after that massacre.
Let us define the events B s (t) and B(t) by In that estimate let us set ε 2 = 0.5 and ε 1 = 1. Taking  where R 4 ,R 5 are defined by (4.40) and (4.41).
To prove (4.35), we set ε 3 = ε 0 /2. Taking into account condition (4.75) for ε 0 we see that if the both events Not A mut,s (t) and Not A sur,s,ε 0 /2 (t) take place, then inequality (4.35) is fulfilled. Thus Finally, taking into account the results of steps 1, 2 and 3 we see that estimate (4.35) holds with the probability Pr t,+ . It completes the proof of (4.35). The second inequality (4.42) can be obtained in the same way.

Remaining part of the proof of Theorem 3.2
We use the same idea that in the proof of Theorem 3.1 but first we establish uniform bounds for the population size and other quantities involved in the proof.
Step 1 Here we estimate the population size. Let us set where q i , q  are defined by relations (4.68)-(4.72). We can find asymptotics of ρ under natural assumptions that p mut → 0 and N popmax → ∞ while all the rest parameters are fixed. Then the leading term in the right hand side of (4. 36 Now repeating the same arguments that in the end of the proof of Theorem 3.1, and taking into account asymptotics (4.77), we obtain the conclusion of Theorem 3.2.

Discussion
In this paper, we proposed a model for fitness landscape learning, which extends earlier work by 7-9 in two ways. First, we use hybrid circuits involving two kinds of variables. The first class of variables are real valued in the interval (0, 1) and can be interpreted as relative levels of phenotypic traits, other variables are Boolean and can be interpreted as genes. Second, we use a threshold scheme of regulation, which is inspired by ideas of the paper by 20. All variables are involved in gene regulation via thresholds.
The work presented here is a major extension of a long term effort to explicitly model the effects of phenotypic buffering in evolution by considering a class of Boolean and mixed Booleancontinuous models in which the phenotype is represented explicitly and the degree of phenotypic buffering can be controlled in various ways. For example, we have demonstrated that the idea of an "evolutionary capacitor" 42,43 can be implemented by explicit control of phenotypic buffering in a hub-and-spokes architecture 23 and that in a more general class of genetic architecture numerical simulations show that an intermediate level of buffering is optimal for evolution in a changing environment 31 .
The results reported here are very promising, since they are consistent with the results of recent experiments by 44 and 45 on heat shock stress. The essential mechanism is that the exploration of the fitness landscape by the genetic network in such a way that future mutations are more likely to be adaptive. We have shown that, at least for some fitness landscapes, rapid evolutionary changes-perhaps instances of the "hopeful monsters" of Goldschmidt 46 -can be created by a combination of random small mutations and epigenetic effects. The main idea is that small mutations pave the way for large epigenetic or genetic changes. The hypothetical mechanism, which we propose, can be outlined as follows (see Figure 4, Figure 5).
This network modifies the thresholds by a simple feedback mechanism. Although we are not aware of clear experimental evidence for the existence of such a mechanism, we nevertheless think that such a mechanism can be connected with regulations via enhancers 47,48 , where enhancer action is described by deep network models based on thermodynamics, and chemical kinetics, and those models contain threshold parameters. Alternative variants involve modifications of weights. To some extent, both mechanisms are mathematically equivalent. However, the regulation via thresholds has an important advantage: it makes phenotypes robust with respect to mutations. In this second version of the paper we included a subsection about regulation and two pictures, which show results of numerical simulations based on the "strong selection weak mutation" (SSWM) algorithm. As is the case for threshold signs, they of course should have different signs. The key question is how an individual obtains information about the fitness landscape. In our model that information is the sign of b j . If the b j is positive the corresponding threshold must be negative, otherwise, it should be negative. Actually, we think that the gene of individuals have that information! According to the main theorems proved in this paper, the fact that a mutated individual survives for a sufficient long period of time gives us that information. The very fact of the existence of the individual carries the most important information.  First we make small random mutations (shown by red vectors), which explore the fitness landscape. If such a mutation is not eliminated from the population, this means that a correct evolution direction is found, and gene regulation system makes a big leap (shown by the green vector) in the direction of that small mutation. Such a two step model can be called clever Goldshmidt leaps. Note that evolution is gradual, and the existence of clusters of almost identical genes involved in the same QTL increases the chances to create a clever Goldschmidt hopeful monster.
Look at the mutations of its genotype and you will know where evolution should go! As is described in Subsection 2.5, we assume that expression of genes involved in the expression of phenotypic traits depends on threshold parameters h j , which take three values: a large negative one, a neutral value close to zero and a large positive one. First the threshold parameter h j is small and thus the phenotypic trait is sensitive with respect to even small mutations. Those mutations play a fundamental role working as scouts exploring environments (see Figure 4). If a mutation occurred and the corresponding mutant has survived within T c  1 generations then according to Theorem 3.1 and Theorem 3.2 these events mean that that mutation increases the fitness that allows the network to estimate the correct direction of evolution. Then gene regulation detects that increase to change the threshold according to simple rules. Namely, if the trait is less expressed in that mutant with respect to wild type parent, the gene regulation system decreases the threshold up to the large negative value. On the contrary, if the trait is strongly expressed in the mutant, the gene regulation system increases the threshold up to the large positive value. This simple regulation control not only sharply reduces the number of mutations needed for adaptation, but also canalizes the phenotype since for large thresholds the trait expression level becomes insensitive with respect to mutations. We suppose that these threshold modifications can be inherited.
So, we propose the mechanism: small mutations serve as scouts finding the way for large epigenetic or genetic changes, which can be performed by gene regulatory system.
The mechanism may also explain the results of 4 on prediction of environmental changes. In fact, let us suppose that environment varies in time. The first, perhaps relatively small, variations can trigger the threshold mechanism described above. As a result, the population will be adapted to the subsequent changes in advance.
Our results show that evolution can proceed rapidly because it reduces the number of mutations required for adaptive change.
The primary limitation of our results is that the representation of the evolving genetic network is limited to the network of gene controlling phenotype, represented here by the Boolean strings s. Other model variables represent the coarse-grained activities of genes. One class is the terminal differentiation genes represented by w ij , and another are the genes or epigenetic factors controlling the thresholds h and their associated learning rules. A more careful consideration of the relationship of these moieties to observable molecular entities is an important objective of future work. At the mathematical level, the key analytical results were obtained in a simplified context that falls short of a realistic level of pleiotropy and thus of the level of NP-hard complexity exhibited by fully pleotropic forms of our model. We believe that our analytical results can be generalized, which we plan to address in future work.

Data availability
All data underlying the results are available as part of the article and no additional source data are required. In the current work the authors study a population genetics model in which fitness is a linear function of a set of phenotypic traits, and where the genotype-to-phenotype map is given by a linear transformation composed with sigmoidal functions. Despite the seeming simplicity of the fitness function, the authors make the case that optimizing fitness is a hard, NP-complete, problem. Under this model, they study the extent to which the fitness landscape (that is, the question of which phenotypic traits contribute positively to fitness, and which contribute negatively) can be inferred from the distribution of these traits in the population after being subject to evolution for a moderate span of time. They then connect this question with that of whether learning algorithms (potentially epigenetic in nature) can help optimize and speed up evolution by learning the fitness landscape.

I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
Toward the goal of inferring the fitness landscape, the authors prove two theorems. Theorem 3.1 concerns a simplified model of an infinite population with no mutations, while Theorem 3.2 concerns a more complex model of a finite population with mutation, stochastic number of descendants and a culling process. In either case, the result is that under certain assumptions, if a mutant genotype is present in the population with high enough frequency after a long enough period of elapsed time, then we can confidently infer that any phenotypic trait differential between the wildtype and the mutant is associated with a higher fitness.
While Theorem 3.1 ignores mutation, even Theorem 3.2 seems at odds with mutation-selection balance. Even a mutant with lower fitness will be present in the population at a frequency that is on the order of the rate of mutation, while the theorem seems to claim that the frequency of such a mutant will be exponentially small with high probability. The authors should resolve this apparent discrepancy. In addition, the main ideas of the proof would be more clearly communicated if the authors would include a treatment of the intermediate model of an infinite population with mutation.
As for the discussion of the way learning can speed up the evolutionary process, this part of the paper remains unclear and underdeveloped. The authors discuss a two-step evolutionary process in which the first step consists of small mutations in order to explore the fitness landscape, and the next step involves changing the thresholds involved in the genotype-to-phenotype map in a way that promotes phenotypic changing the thresholds involved in the genotype-to-phenotype map in a way that promotes phenotypic traits associated with higher fitness. While this idea is interesting and worth exploring, a few issues arise. A conceptual issue remaining to be addressed is whether the threshold h is part of the genotype and what mechanism are needed to alter its value. According to the authors, the h 's can be modified genetically or epigenetically. If epigenetically, it is not apparent what are the environmental cues that will lead to such learning let alone the actual mechanism of modifying them. If genetically, it is similarly unclear in what way the learning of the fitness landscape is being stored, if at all, in the genotype, and what is the connection to the Theorems of chapter 3. The theorems in chapter 3 rely on observing the frequency of a genotype in the population, but such information is not stored in individual genotypes. Additionally, if the thresholds are understood to be variable and subject to selection, then the fitness-maximization problem in fact becomes easy (in contrast to the prior analysis of it as NP-complete), unless we impose restrictions on the range of the threshold. Indeed, one can set the thresholds at positive or negative infinity depending on whether the corresponding trait is positive or negative in order to effectively keep the trait on or off regardless of genotype.
Finally, as this manuscript addresses the relation between learning and the rate of evolution, it would benefit from including a reference to one of the most relevant and intuitive articles written in the late 80's by Geoffrey E. Hinton & Steven J. Nowlan, "How Learning Can Guide Evolution" , 1, Complex Systems 495-502 . In it, Hinton and Nowlan showed how learning alters the shape of the fitness landscape and thereby provides easier evolutionary paths towards sets of co-adapted alleles. Hinton and Nowlan demonstrated that this effect allows learning individuals to evolve faster than non-learners. Though the learning model presented by Hinton and Nowlan operates at "somatic" timescale, the analogy to mutations at the evolutionary timescale can be drawn.
To conclude, though it can be improved by the above suggestions, this article touches upon very interesting and important issues in the field of evolutionary biology which are still only lightly investigated, and highlight what might be a fruitful path towards better understanding.
No competing interests were disclosed. would in this case be nonnegative) then I do not understand the thought example. Increasing the threshold would reduce the expression level that was selected to go up in the first place! These points seem to me crucial, so their clarification is badly needed. I would appreciate a few numerical examples along with at least hints to answers to the more conceptual questions on the dynamics. Also, according to the cited Nagylaki dynamics the population evolves almost as if it were in linkage equilibrium -which cannot hold for the asexual population considered by the authors. There are also some minor issues in the paper: in Eq. 2.4 "s" should be replaced by a sigma, and there are also examples where plurals and singulars in the same sentence do not match.
Your article is published within days, with no editorial bias You can publish traditional articles, null/negative results, case reports, data notes and more The peer review process is transparent and collaborative Your article is indexed in PubMed after passing peer review Dedicated customer support at every stage For pre-submission enquiries, contact research@f1000.com