Adaptation, fitness landscape learning and fast evolution

2.1 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 between Boolean hypercube and genotypes is considered for example in 12.

2.2 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₁, . . . , f_Nb) 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₁, i₂, . . ., 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 S_g = {1, . . . N_g}

We denote the total number of genes involved in regulation of all f_j by N_r, where $N_r=\sum n_j\le K{N}_b\mathrm{.}$

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

where j = 1, . . . , N_b. Here σ(z) is a sigmoidal function of real z such that

and w_ij, h_j are some coefficients (their meaning will be explained below). As an example, we can take σ(S) = (1 + exp(−bS))⁻¹, 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²⁰. 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 $n_j^{\pm }$, respectively.

Then f_j ≈ 1 if $n_j^{+}-n_j^{-}\gg h_j\mathrm{.}$

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.

2.3 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²⁴. 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²⁶, 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²⁷. In evolution computations (EC) some artificial fitness models were used, for example OneMax and Leading Ones to test evolution algorithms, see for example¹⁵.

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

where K_F is a positive constant and b = (b₁,..., b_Nb) 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 (2.7) can describe gene epistatic effects via dependence of f_j on s if f_j are nonlinear in s.

The expression (2.7) can serve as a rough approximation of the fitness function in the case of insects such as grasshoppers or fruit flies. In fact, important factors, which determine insect survival, are thermoregulation, mimicry and camouflage levels^18,22,28. All those factors depend on colour pigmentation pattern. Blackwhite pigmentation patterns can be roughly described by vectors f = (f₁, f₂,..., f_Nb), where f_j ≈ 1 and f_j ≈ 0 mean that the cell j is black, or white, respectively, Then thermoregulation depends on ${\displaystyle {\sum }_j{f}_j}.$ The mimicry level can be approximately defined by expression ${\displaystyle {\sum }_j\left|f_j-f_j^{*}\right|}$ , 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¹⁸. Second, the topology of the fitness landscapes was studied in 22 by field experiments in the case of insect mimicry. Main conclusions are as follows. A number of studies of fitness landscapes in natural populations have demonstrated low fitness of intermediate phenotypes, i.e., existence of valley in the fitness landscape. It is found²² that natural selection promotes genetic architecture preventing the expression of intermediate phenotypes. Close fitness peaks are separated by ridges, favouring colour pattern switches and allowing drift from local peaks.

In Section 2.6 we will show that the fitness model defined by (2.4) and (2.7) have those topological properties.

2.4 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:

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

where $\overline{F}$(t) is the average fitness of the population at the moment t defined by

where S(t) is the set of genotypes existing in the population at time t (the genetic pool) and X(s, t) = N(s, t)/N_pop(t) is the frequency of the genotype s. Here N(s, t) denotes the number of the population members with the genotype s at the step t.

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 $\overline{F}$(t) defined by (2.10) satisfies Fisher’s theorem, so that this function increases at each time step t: $\overline{F}$(t + 1) ≥ $\overline{F}$(t).

2.5 Gene regulation network

In this section, we follow ideas of the classic paper²⁹: 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 $h_i^{(-)}$, $h_i^{(+)}$ and $h_i^{(0)}$ such that

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: if S_i > 0 is large enough and h_i is small, then h_i = −D_h; if S_i < 0 and $\left|S_i\right|$ is large enough and h_i is small, then h_i = D_h; otherwise, we do not change h_i.

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 $\left|S_i\right|$.

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 $\left|h_i\right|$ are non-sensitive with respect to mutations. The effect produced by this robustness is shown on Figure 1.

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.

This regulation is even more effective, if we consider the modular evolution, following the recent paper³². 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 $\left|h_i\right|$, 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.

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

Figure 3. Frequency distributions of degree of gene pleiotropy for model (2.4) with the parameters N_g = 4000, β = 4, h = 0, N_b = 3000.

2.6 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. 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³³ 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

As a consequence of the second assumption, F attains its maximum for f₁ = 1, f₂ = 1,..., f_Nb = 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.

K-SAT problem. Let us consider the set V_n = {s₁,..., s_n} of Boolean variables s_i ∈ {0, 1} and a set 𝒞_m of m clauses. The clauses C_j are disjunctions (logical ORs) involving K literals z_i1, z_i2,..., z_ik, where each z_i is either s_i or the negation $\overline{s}$_i of s_i. The problem is to test whether one can satisfy all of the clauses by an assignment of Boolean variables.

Cook and Levin^34,35 have shown that the K-SAT problem is NP-complete 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 N_b < α_c(K)N_g, where α_c(K) ≈ 2^K log 2 for large K.

The set 𝒞_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 𝒞_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)³⁹. 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.

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 $S_j={\displaystyle {\sum }_{i=1}^{N_g}{w}_{ji}{s}_i\mathrm{.}}$ We set m = N_b and n = N_g. Under this choice of h_j, the terms σ(S_j) can be represented as disjunctions of literals z_j. Each literal z_j equals either s_j or $\overline{s}$_j, where $\overline{s}$_j denotes negation of s_j. To 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⁴¹, 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 𝒞_K, mentioned above, outline the properties of really observed fitness landscapes²²: 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.

3.1 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 ∊ I_– expression of the trait decreases the fitness.

Let s and $\overline{s}$ be two genotypes. Then we denote by Diff(s, $\overline{s}$) the set of positions i such that s_i ≠ $\overline{s}$_i:

The set Diff(s, $\overline{s}$) indicates which genes in s should be flipped in order to obtain $\overline{s}$.

We formulate two theorems on fitness landscape learning. First we consider the case of infinitely large populations.

Theorem 3.1. Suppose that the evolution of the genotype frequencies X(s, t) is determined by equations (2.9) and (2.10). Moreover, assume that

I for all t ∊ [T₁, T₁ + T_c], where T₁, T_c > 0 are integers, the population contains two genotypes s and $\overline{s}$ such that the frequencies X(s, t) and X($\overline{s}$, t) satisfy

II we have

for some j. In other words, the genes s_i such that s_i ≠ $\overline{s}$_i are involved in a single regulation set K_j; and finally,

III Let

and

Then, if

we have j ∊ I₊. If f_j(s) > f_j($\overline{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 $\overline{s}$ are contained in a single regulatory set K_j.

The assertion of Theorem 3.1 is not valid if the set Diff(s, $\overline{s}$) belongs to a union of different regulation sets K_j , j = j₁, . . . , 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 $\overline{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₁

Consider inequality (3.11) for T = T₁ + T_c. Let us note that in the relation Q(T₁) = X(s, T₁)/X($\overline{s}$, T₁) the numerator is p₀ whereas the denominator ≤ 1. Thus, Q(T₁) ≥ p₀. The same arguments show that Q(T₁ + T_c ) ≤ 1/p₁. Therefore, by (3.11) one obtains that

This inequality leads to a contradiction for T_c satisfying (3.6), thus completing the proof.

3.2 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

where c_F, C_F > 0 are constants independent of t. For example, if

then c_F = K_F exp(–γ) and C_F = K_F exp(γ) and (3.13) holds if K_F > exp(γ).

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.

Let

Note that for simplicity in the next Theorem 3.2 we consider point mutations (bit flipping) only. The model used here cannot represent mutations of arbitrarily small effect, but it can include insertions or deletions. In contrast to Theorem 3.2, Theorem 3.1 is valid for all kinds of mutations.

Then we have

Theorem 3.2. Consider the population dynamics defined by model 1-3 in Subsection 2.4. Assume conditions (3.14) and M hold, and assumptions (3.3), (3.4), (3.5), (3.7) of Theorem 3.1 are satisfied. Suppose

Then if j ∊ I_– the inequality

is fulfilled with the probability Pr_v such that

where for large N_popmax and p_mut → 0

Interpretation of Theorem 3.2

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 H_– that j ∊ I_– against the hypothesis H₊ that j ∊ I₊. Suppose that H_– is true. Let V be the event that the frequency X($\overline{s}$, t) of the genotype $\overline{s}$ in the population is larger than p₁ 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 H_– should be rejected. We will refer T_c as the checking time.

Rare mutants. In this Theorem we assume that the frequencies p₀ and p₁ 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 $\overline{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.

are satisfied, where F(s) is a fitness of genotype s, $\overline{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.

4.1 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^{\prime }$ (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^{\prime }$(x, t) is

where EX denotes the expected value of X. By $\stackrel{-}{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

Moreover, we use two elementary estimates. Let 𝒜 be an event in stochastic population dynamics. We denote by Not𝒜 the negation (complement) of 𝒜 and by Pr(𝒜|ℬ) the conditional probability of 𝒜 under the condition ℬ. For events 𝒜 , ℬ₁, . . . , ℬ_n we have

For two events 𝒜, ℬ one has

Lemma 4.1. Let X_i be independent random quantities, where i = 1, . . . , n. Let each X_i be distributed according to the Poisson law with the average EX_i = μ_i. Let us denote

Then for all δ > 0

where

Similarly,

Proof. Note that for any λ > 0

Since X_j are independent quantities, we have

The straight forward computation shows that

Therefore, due to the Markov inequality (4.3) and estimate (4.8) one has

where

We minimize f with respect to λ and obtain (4.6). To derive (4.7), we use

and repeat the same arguments. The Lemma is proved.

Lemma 4.2. Let X_i be independent random quantities, where i = 1, . . . , n such that X_i ∊ {0, 1} and EX_i = p. Then

where

Proof. Note that for any λ > 0

Since X_j are independent quantities, we have

Note that E exp(−λX_j) = p exp(−λ) + 1 − p. Let

We take λ = ln 2 and find that G(ln 2, p) = −g(p). Now by using the Markov inequality (4.3) and estimate (4.12) one obtains (4.10). The Lemma is proved.

We also use the following Chernoff-Hoeffding theorem. Let X_i be i.i.d. quantities such that X_i ∊ {0, 1} and EX_i = p, where i = 1, . . . , n. Then for $X={\displaystyle {\sum }_{j=1}^n{X}_j}$ one has

where D(x||y) is the Kullback-Leibler divergence

Moreover, we will use the Hoeffding Theorem: if i.i.d. quantities X_i ∊ [0, 1] with the probability 1 then

4.2 Main lemmas

First we estimate the population size fluctuations.

Lemma 4.3. Let $\overline{N}$(t) be the number of all progeny, born in the population at the moment t before the massacre, and ε₁ > 0 be a small number. Then

with probability

where

Proof. Let $n^{\prime }$ (x, t) denote the number of progeny produced by the individual x before the massacre at the t-th evolution step. The number $\overline{N}$(t) is the sum

of the mutually independent random quantities. According to (4.2), the average E$N^{\prime }$ (x, t) is F(s_g (x)). Therefore,

We set

and use the Lemma 4.1 that gives us (4.17).

Lemma 4.4. Let ε₂ ∊ (0, 1) be fixed and condition (3.13) be fulfilled. Assume, moreover, that

where

and c_F > 1 is defined by (3.13). Let us define the event 𝒟_ε2 (t) by

Then one has

where

and

Proof. Let ξ(x) be random quantities defined as follows: ξ(x) = 1 if the individual x is survived as a result of massacre (see point 3 of our model from Subsection 2.4), and ξ(x) = 0 otherwise. Let $X^{\prime }$ (t) be the set of progeny produced by all individuals from the population. Then the number N_sur(t) = N_pop(t + 1) of finally survived progeny can be computed as follows:

Note that |$X^{\prime }$(t)| = $\overline{N}$(t). Moreover, Eξ(x ) = N_popmax/$\overline{N}$(t) for $\overline{N}$(t) ≥ N_popmax. Therefore, if $\overline{N}$(t) ≥ N_popmax then

Let us define the event

where the interval J_ε(t) is defined by (4.16) and $\tilde{\epsilon }$ is defined by (4.27). By (4.4) we have

Now we apply the Hoeffding inequality (4.15). For each ε₂ > 0 we obtain

If ℬ (t) takes place, then $\overline{N}$(t) ≥ N_popmax and consequently

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 $\overline{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.

Lemma 4.5. Let ε₀ be a positive number satisfying (4.75) and

Then one has

with the probability Pr_s,t,+ such that

where

Similarly,

with the probability Pr_s,t,− such that

Proof. Step 1, estimates of fluctuations. First let us estimate the fluctuations of the number $\overline{N}(s,t)$. For each ε₂ > 0 let us define the event

By Lemma 4.1 one has

Note that

As a result, by (4.46) we obtain

Step 2. Here we estimate the number of progeny that survived as a result of the massacre procedure (point 3 of the population dynamics model, see subsection 2.9). Let X' (s, t) be the set of progeny produced by individuals with the genotype s. Then the number N_sur(s, t) of survived progeny x for individuals x belonging to the set Z' (s, t) is

where ξ(x) are defined in the proof of the previous Lemma. For ε₃ > 0 we consider the event

Let us estimate the probability Pr(𝒜_sur,s(t)). According to the Hoeffding Theorem (4.15)

Note that ξ(x) and ξ(y) are independent quantities for different x and y, thus under the condition $\overline{N}(t)$ > N_popmax

therefore,

Let us define the events ℬ_s(t) and ℬ(t) by

Then using (4.4) one has

We observe that under conditions ℬ_s(t) and ℬ(t)

In that estimate let us set ε₂ = 0.5 and ε₁ = 1. Taking into account that E$\overline{N}(t)$ = $\overline{F}(t)$N_pop(t) < 2C_FN_popmax, we have that

where

Moreover, according to (4.47)

and due to (4.17)

where R₁,R₂ are defined by (4.37) and (4.38). Finally,

Step 3, estimate of the number of mutants.

Let us estimate how many individuals with genotypes s can mutate. The probability of mutation is p_mut. Let N_mut(s, t) be the number of such mutants. Let us define the event 𝒜_mut,s(t) by

Since the random quantity N_mut(s, t) is subject to the Bernoulli law, we can apply the Chernoff-Hoeffding inequality (4.13). Then we obtain that

where, according to definition (4.14) of D(x||y), one has

and g is defined by (4.11).

Using (4.4) one has

As a result, by Lemma 4.3 one finds

where R₄,R₅ are defined by (4.40) and (4.41).

To prove (4.35), we set ε₃ = ε₀/2. Taking into account condition (4.75) for ε₀ we see that if the both events Not 𝒜_mut,s(t) and Not 𝒜_sur,s,ε0/2(t) take place, then inequality (4.35) is fulfilled. Thus

where R_i are defined by (4.37)–(4.41).

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.

4.3 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

in Lemma 4.4. Let us consider the events 𝒟_ε2(t) defined by (4.24) in Lemma 4.4. If the events 𝒟_ε2(t) take place for all t ∈ [T₁, T₁ + T_c] and N_pop(0) = N_popmax, we have that