Populations of genetic circuits are unable to find the fittest solution in a multilevel genotype-phenotype map

The evolution of gene regulatory networks (GRNs) is of great relevance for both evolutionary and synthetic biology. Understanding the relationship between GRN structure and its function can allow us to understand the selective pressures that have shaped a given circuit. This is especially relevant when considering spatiotemporal expression patterns, where GRN models have been shown to be extremely robust and evolvable. However, previous models that studied GRN evolution did not include the evolution of protein and genetic elements that underlie GRN architecture. Here we use toyLIFE, a multilevel genotype-phenotype map, to show that not all GRNs are equally likely in genotype space and that evolution is biased to find the most common GRNs. toyLIFE rules create Boolean GRNs that, embedded in a one-dimensional tissue, develop a variety of spatiotemporal gene expression patterns. Populations of toyLIFE organisms choose the most common GRN out of a set of equally fit alternatives and, most importantly, fail to find a target pattern when it is very rare in genotype space. Indeed, we show that the probability of finding the fittest phenotype increases dramatically with its abundance in genotype space. This phenotypic bias represents a mechanism that can prevent the fixation in the population of the fittest phenotype, one that is inherent to the structure of genotype space and the genotype-phenotype map.


58
These Boolean GRNs are built on top of a simple model of cellular biology, t o yLIFE [25,26]. t o yLIFE is enough to steer evolving populations towards more abundant GRNs, thus introducing an additional element when trying to explain GRN evolution, one that is not related to function or structure. Fur-70 thermore, we also show that this phenomenon can result in the inability of the evolutionary search to 71 find some regulatory patterns, even when they are fitter than every other.
In a Boolean GRN, a gene can either be ON (expressed) or OFF (not expressed). The expression state (c) A double-positive feedback loop with self-activation loops leads to both proteins A and B (grey) being expressed in the tissue in a stable way, and expanding through the tissue. (d) A double-positive feedback loop without self-activation leads to an alternating pattern where the tissue expresses first protein A (orange), then protein B (blue), and so on. Notice how the speed with which the pattern extends throughout the tissue is half the speed of patterns in a and b. This is because only protein A is allowed to propagate to the neighbouring cells (see Figure 2), so that the pattern can only extend when protein A is expressed.
negative feedback loop with self-activation, (b) the same as before but with gene a having constitutive expect evolutionary dynamics to choose among these sixteen GRNs with equal probability, everything 124 else being equal. This is certainly what almost every mathematical model of phenotypic evolution 125 (including previous models of GRN evolution) would predict. 126 We performed Wright-Fisher evolutionary simulations with t o yLIFE organisms in a strong selection, 127 weak mutation regime (Methods), and selected the pattern in Figure 3b as the evolutionary target -i.e.

128
we assigned maximal fitness to it, and every other pattern became less fit as it differed more from the 129 target (see Methods for the complete definition of the fitness function). We found that, after 100, 000 130 mutations, 93% of simulations ended up finding one particular GRN among all sixteen (GRN XI in Figure   131 S9, see below), and the network in Figure 3b (GRN V) does not appear as the endpoint of evolutionary 132 dynamics in any of the 1, 000 simulations. In order to understand this somewhat unexpected result, we 133 now discuss how Boolean GRNs are obtained from t o yLIFE genotypes.

134
Regulation in t o yLIFE 135 We will introduce gene regulation in t o yLIFE through an example (for an in-depth discussion of t o yLIFE's genotype given each input, i.e. its truth table (Figure 4d). When no protein is present, the polymerase (which is always present in the cell) will activate gene b and the output will consist of protein B. as a result of protein products propagating from one cell to the next. With this information we can 155 unequivocally compute each genotype's corresponding cellular automaton.

156
It is worth noting that in the process of defining these phenotypic levels we have already introduced 157 a lot of degeneracy. For instance, there are 2 40 ≈ 10 12 genotypes with two genes, but they only give 158 rise to 1, 472 different GRNs, which in turn generate only 453 different cellular automata -an average 159 of ≈ 2 × 10 9 genotypes per cellular automata. Not all GRNs are equally probable in genotype space, 160 however: the distribution of abundances of GRNs follows a log-normal distribution (Figure 4e), which 161 has been observed in many other genotype-phenotype models and has been shown to be universal under 162 some very general assumptions [16,[26][27][28]. The most abundant GRN is mapped by more than 500 163 billion genotypes, while the rarest one is only mapped by 8 genotypes. This phenomenon has been 164 called phenotypic bias [15,16], and it is also observed in the distribution of abundances of cellular automata (Supplementary Figure S10a). As a consequence, the sixteen GRNs that generate the pattern   Figure S11). In comparison, the least abundant Boolean GRN among these sixteen (GRN VIII) is Pattern 113 is rarely found in our evolutionary simulations. b) Pattern 109 is similar to 113, but it is generated by 1.64 × 10 8 genotypes -about 10 5 times commoner. As a result, it appears as the endpoint of our simulations 84% of the times. c) Pattern 170 also appears as the endpoint of the simulations 8% of the times, even though it is not very similar to pattern 113. This is due to its high abundance in phenotype space: 1.36 × 10 8 genotypes are mapped to it. (d) This phenomenon is not restricted to pattern 113. The probability of finding a target pattern (p) goes to zero as the logarithm of pattern abundance (S) decreases. Line: p = (1 + (430767/S) 1/2 ) −1 , R 2 = 0.58. (e) Even when simulations do find the fittest pattern, the time to reach it (T ) increases as pattern abundance decreases. Line: log 10 T = 4.35 − 0.05(log 10 S − 2.93) 2 , R 2 = 0.68.
it appears as the evolutionary endpoint only in 3% of the 1, 000 simulations. Instead, the pattern that 191 appears in most of the simulations is pattern 109 (Figure 5b), which has a fitness of 0.991 relative to 192 that of pattern 113, and is mapped by 1.6 × 10 8 genotypes. There are 33, 280 mutational paths between The fitness function for our evolutionary simulations is calculated as follows: each pattern is a string in 262 base four of length L = 31 · 100. For every evolutionary scenario, we choose one particular pattern p T 263 as the target value, and assign fitness 1 to it. Then we compute the Hamming distance D of a pattern 264 p to the target as where d i, j is Kronecker's delta, which is equal to 1 if i = j and 0 otherwise, and p(i) is the i-th letter in 266 the string p. Fitness is then calculated as Evolutionary simulations 268 We assume a strong selection, weak mutation scenario. In this regime, Wright-Fisher dynamics are 269 reduced to a continuous-time random walk in genotype space. We only consider point mutations, which 270 arise in the population at constant rate µ, and the fixation rate of a new mutation is given by where f is the fitness of the current phenotype relative to that of the mutant, and N is population where n = |x| and N w (x) is the number of words in the dictionary created by the Lempel-Ziv algorithm  sequence of length 16, we fold it into every SAW and compute its folding energy, following the HP model. For instance, we fold the sequence PHPPPPPPPPPHHHHP into one of the SAWs and compute its folding energy (right). There are two HH contacts, five HP contacts and two PP contacts -we only take into account contacts between non-adjacent toyAminoacids. Summing all this contacts with their corresponding energies, we obtain a folding energy of −11.5. Repeating this process for every SAW, we obtain the minimum free structure.
reflections-, there are only 38 SAWs on that lattice (Supplementary Figure S2).

417
The energy of a fold is the sum of all pairwise interaction energies between toyA that are not 418 contiguous along the sequence. Pairwise interaction energies are E HH = −2, E HP = −0.3 and E PP = 0, 419 following the conditions set in [40] that E PP > E HP > E HH (Supplementary Figure S2). toyProteins 420 are identified by their folding energy and their perimeter. If there is more than one fold with the 421 same minimum energy, we select the one with fewer H toyAminoacids in the perimeter. If still there is 422 more than one fold fulfilling both conditions, we discard that protein by assuming that it is intrinsically 423 disordered and thus non-functional. Note, however, that sometimes different folds yield the same folding 424 energy and the same perimeter. In those cases, we do not discard the resulting toyProtein. In those cases we select the sites that yield the most centered interaction (Supplementary Figure S4b).

490
If ambiguity persists, no bond is formed. Also, no more than one toyProtein / toyDimer is allowed to 491 bind to the same toyMetabolite, even if its length would permit it. toyProteins / toyDimers bound to 492 toyMetabolites cannot bind to promoters. and two copies of P2 will form two copies of dimer P1-P2, and one copy of P1 will remain free.

516
Expression of toyGenes occurs through the interaction with the toyPolymerase, which is a special kind 517 of toyProtein (see Supplementary Figure S1). The toyPolymerase only has one interacting side (with 518 sequence PHPH) and its folding energy is fixed to value −11.0: it is more stable than more than half 519 the toyProteins. It is always present in the system. The toyPolymerase binds to promoters or to the 520 right side of a toyProtein / toyDimer already bound to a promoter. When the toyPolymerase binds to a 521 promoter, translation is directly activated and the corresponding toyGene is expressed (Supplementary 522 Figure S5a). However, a more stable (lower energy) binding of a toyProtein or toyDimer to a promoter when the toyPolymerase binds to its promoter region. The sequence of Ps and Hs of the toyProtein will be exactly the same as that of the toyGene coding region. (b) If a toyProtein binds to the promoter region of a toyGene with a lower energy than the toyPolymerase does, it will displace the latter, and the toyGene will not be expressed. This toyProtein acts as an inhibitor. (c) The toyPolymerase does not bind to every promoter region. Thus, not all toyGenes are expressed constitutively. However, some toyProteins will be able to bind to these promoter regions. If, once bound to the promoter, they bind to the toyPolymerase with their rightmost side, the toyGene will be expressed, and these toyProteins act as activators. (d) More complex interactions -involving more elements-appear. For example, a toyProtein that forms a toyDimer with an inhibitor -preventing it from binding to the promoter-will effectively activate the expression of the toyGene. However, it does neither interact with the promoter region nor with the toyPolymerase, and its function is carried out only when the inhibitor is present.
We call this kind of toyProteins conditional activators. (e) Two toyProteins can bind together to form a toyDimer that inhibits the expression of a certain toyGene. As they need each other to perform this function, we call them conditional inhibitors. As the number of genes increases, this kind of complex relationships can become very intricate.
act as activators (Supplementary Figure S5c). This process finds a counterpart in toyProteins that bind to promoter regions more stably than the toyPolymerase does, and therefore prevent gene expression -533 this happens if E int(PROT) + E PROT < E int(POLY) + E POLY . They are acting as inhibitors (Supplementary 534 Figure S5b). There are two additional functions that could not be foreseen and involve a larger number 535 of molecules. A toyProtein that forms a toyDimer with an inhibitor -preventing its binding to the When a toyDimer is bound to a toyMetabolite, another toyProtein can interact with this complex and 546 break it. This reaction will take place if the toyProtein can bind to one of the subunits of the toyDimer 547 and the resulting complex has less total energy than the toyDimer. As with the rest of interactions, 548 the catabolic reaction will only take place if this binding is unambiguous. As a result of this reaction, 549 the toyDimer will be broken in two: one of the pieces will be bound to the toyProtein (forming a new 550 toyDimer), and the other one will remain free. The toyMetabolite will break accordingly: the part of 551 it that was bound to the first subunit will stay with it, and the other part will stay with the second 552 subunit. Note that the toyMetabolite need not be broken symmetrically: this will depend on how the a new toyDimer energetically more stable than the old one, the two toyProteins will unbind and break the toyMetabolite up into two pieces. We say that the toyMetabolite has been catabolised. .216 x 10 9 II) 10 6 III) 1.652 x 10 9 IV) 3 x 10 6 V) 9 x 10 6 VI) 1.71 x 10 8 VII) 3.0 x 10 7 VIII) 0.2 x 10 6 IX) 3.536 x 10 9 X) 6 x 10 6 XI) 2.017 x 10 11 XII) 1.9 x 10 7 XIII) 1.42 x 10 8 XIV) 10 6 XV) 4.91 x 10 8 XVI) 1.293 x 10 9 Supplementary Figure  patterns can also be fitted by a log-normal distribution, although the fit is rather noisy (R 2 = 0.41), given that we only have 172 patterns to fit.  find every GRN with equal probability. On the contrary, those GRNs that are more abundant in genotype space appear more frequently as an endpoint of our simulations, in agreement with Refs. [15,30]. In fact, the fraction of times a given GRN is the endpoint of the simulations is almost exactly its abundance in genotype space relative to that of all sixteen GRNs. Linear fit is approximately y = x (R 2 ≈ 1.0).