SeRANN model

Formally, a SeRANN is an artificial neural network that approximates a function , such that

The inputs are X, a grayscale image with m = 28 rows and n = 28 columns of pixels, and g, a genotype bit-string with k = 100 bits (i.e. zero or one values). The outputs are , a probability vector of size L = 10, where each entry in this vector represents the probability that the image X is of the i-th class (); and g′, a bit-string genotype, also with k bits. See Fig 1 for an illustration.

The objectives of the SeRANN are (i) to accurately classify X to the correct class, by giving the highest probability for the correct class i, and (ii) to faithfully replicate g, such that g′ is as similar as possible to g. A SeRANN is trained using stochastic gradient descent with backpropagation, a standard practice in deep learning [22] with the Adam optimizer [23]. The network’s (trainable) parameters are randomly initialized before training using Glorot initialization [24]. This initialization scheme is commonly used to prevent exploding and/or vanishing gradients and to accelerate training convergence. We define two loss functions. First, ℓ_X is the categorical cross entropy between the ground-truth vector y and the predicted probability vector ,

Second, ℓ_g is the mean squared error between the input genotype g and the output genotype g’, where g_j is the j-th bit of genotype g. The loss function ℓ used for training SeRANNs is a weighted average of the above loss functions, (1) where the weight 0≤α≤1 is determined by a SeRANN source-code variable, loss_weight. Hence, mutation and selection can change the value of this variable, which effectively controls how much of the training effort is dedicated to each of the two tasks.

Source code representation

We implement SeRANN using the Python programming language [17] and the Keras deep-learning library [19]. This library provides a simple API for defining, training, and evaluating artificial neural networks. In general, each line describes a layer in the network, including the layer type, its hyperparameters and the previous layer it is added upon.

The full SeRANN source code includes the loading of necessary packages (import statements); declaring constant experiment parameters (classification image shape, genotype length); the SeRANN layers definition; and the initialization of the network model object. Of these, only the SeRANN layers definition varies between different SeRANN individuals. Therefore, the SeRANN genotype only encodes this part of the source code, see lines 9–14 in Fig B in S1 Text, which shows a full example source code of a SeRANN implementation.

Genetic code and ribosomal autoencoder

In biology, the genetic code is used to translate (i.e. decode) information encoded in the genetic material (DNA or RNA molecules) to proteins, which, through interactions with the biotic and abiotic environment, give rise to the phenotype. This translation is done by the ribosome. Likewise, we use a genetic code to translate a bit-string genotype g to a Python source-code string (Fig B in S1 Text). This source code is then interpreted by the Python interpreter to an instance of a neural network model implemented in Keras, and trained by supervised learning, gradient descent, and backpropagation, to produce the”adult” phenotype.

The viable SeRANN phenotypes space is constrained: SeRANN source code must be exe- cutable by Python interpreter. This means that the genotypes space topology and the genotype- phenotype encoding must be designed with a special care. Otherwise, the vast majority of mutants (i.e. errors in replicated genotypes) would be doomed to immediate “death”, as their translation will suffer from invalid syntax and fail to execute. The fitness landscape induced by such flawed encodings would be sparse, or”perforated”—many potential solutions will be virtually unreachable and the efficiency of the evolutionary process will reduce [25]; see also “holey” adaptive landscapes [26].

To translate genotypes to source codes, we introduce the ribosomal autoencoder (RiboAE). The RiboAE is an entirely separate neural network from SeRANN. We have trained it to encode Python source code to a genotype bit-string, and faithfully decode it back to source code, using unsupervised learning, which only requires examples of source codes, without examples of genotypes. After training the RiboAE, it was kept constant, and did not change during SeRANN evolution.

SeRANN source-code generator

To generate a significant number of synthetic SeRANN source-code examples, sufficient for training the RiboAE, we implemented a heuristic SeRANN source code generator. This generator is a stochastic state machine (i.e., Markov chain) that defines the probability of transition between a limited set of manually selected ANN layers (see Table D in S1 Text). For example, given that the current layer is a convolutional layer, the next layer will be a pooling layer with probability 0.7, another convolutional layer with probability 0.2, or a fully-connected layer with probability 0.1.

The hyper-parameters of each network layer were drawn from a discrete normal distribution around reasonable expected values. For example, we sampled the number of neurons in a fully-connected layer from a normal distribution, N(μ = 64,σ = 15), clipped the values to lay between 8 and 128, and rounded them to the nearest integer.

To ensure the validity of the generated SeRANN source codes, we enforced a general structure: a fixed-size input layer, a merge layer (to join both classification and replication input branches) and a fixed-size output layer. Thus, we allowed SeRANN source-code examples to diverge by the number of hidden layers, and their types and hyper-parameters, while keeping them compatible with the replication and classification tasks. We used the SeRANN generator to generate 1,000,000 examples of SeRANN source codes.

RiboAE architecture and training

The RiboAE takes a neural network source code as input and is trained to output the exact same source code (Fig A in S1 Text). The d-dimensional input and output of the RiboAE is a SeRANN Python source-code string (Fig B in S1 Text), segmented into a sequence of atomic language tokens. These tokens include variable names and literals (digits, brackets, etc.). The sequence of tokens is padded to exactly d tokens using a special padding token. Tokens are mapped to integers according to a predefined vocabulary (see complete vocabulary in Table A in S1 Text).

The RiboAE is composed of two components: the RiboEncoder and the RiboDecoder. The RiboEncoder consists of an embedding layer, three convolutional layers, and a fully connected layer. It takes as input a d-dimensional sequence that represents the source-code string as d tokens, and encodes it to a binary k-dimensional genotype vector g, where k ≪ d. This is parallel to the molecular process of reverse translation (note that reverse translation and reverse transcription are different processes). The RiboDecoder consists of a single convolutional layer and a fully connected layer. It takes a binary k-dimensional genotype vector g as input and reconstructs the original d-dimensional source-code input, by estimating the categorical probabilities for each of the d tokens. This is parallel to the molecular process of translation.

In practice, the output of the RiboEncoder is a k-dimensional continuous variable, ϕ, where each of its values (ϕ₁, ϕ₂,…,ϕ_k) is in the interval [0,1]. During training, ϕ is not provided directly to the RiboDecoder, but rather interpreted as parameters of a k-dimensional Bernoulli distribution. A random sample from this distribution is used as input to the RiboDecoder.

By training it to overcome minor replication errors, this random sampling procedure allows the RiboDecoder to become robust to mutations. It also reduces the number of lethal “holes” in the genotype space, as reflected by the distribution of fitness effects of new mutations. The fraction of lethal mutations in our main experiment is 10%-20% lower than the one estimated for viruses and yeast. Because the RiboAE is trained using backpropagation and gradient descent [22] with the Adam optimizer [23], the aforementioned random sampling poses two problems: backpropagation through random nodes, and gradient descent with discrete variables, which are not differentiable. We solve the first problem using the reparameterization trick [27]: we transform a sample from a uniform distribution rather than directly sampling from a Bernoulli distribution. Thus, the RiboAE reconstruction error can backpropagate through the RiboDecoder and ϕ to the RiboEncoder.

The second problem arises because we would like our genotypes to be consistent with natural genotypes, which are discrete and consist of combinations of four nucleotides, A, G, C, and T. Here, we use binary genetics (0 and 1), although other genetic systems such as quaternary genetics (A, G, C, and T) are possible in future work. We solve the problem of non-differentiability using the Concrete distribution, a “relaxed” version of the Bernoulli distribution [28]. For each output of the RiboEncoder, ϕ_i, we sample a value from a continuous uniform distribution, u_i ∼ U(0, 1), and use the reparameterization trick, followed by a sigmoid transformation (i.e. logistic function) to produce the genotype value g_i,

The resulting vector, g, consists of continuous values which are mostly very close to 0 and 1. During inference, differentiability is no longer required. Instead of the discrete relaxation, the values of g_i are rounded to 0 or 1,

The RiboDecoder input is the vector g. Its output consists of d vectors of length r, , where r is the number of language tokens in the vocabulary, such that v_i = (v_i,1,…, v_i,r). The softmax transformation is applied to each vector v_i, and the result is interpreted as a parameterization of a categorical distribution over the r possible language tokens (see Table A in S1 Text for the mapping of tokens to integers). Recalling that s is the input token sequence of the RiboAE, and s_i is the integer mapping of the i-th token in s, the probability of s_i given v_i is (2)

The RiboAE is inspired by the variational autoencoder (VAE) of Kingma and Welling [27]. In both cases, the autoencoder is trained to minimize the reconstruction error between the input and output. However, in contrast to the RiboAE, the VAE is simultaneously trained to minimize the Kullback–Leibler divergence between the distribution defined by the encoder output and a prior isotropic distribution (usually a standard normal distribution). Stochastic replication has a biological interpretation: each cell division is equivalent to sampling a genotype from some distribution defined by the original genotype. However, prior constraints on this distribution, as expressed by its divergence from some parametric distribution, have no clear biological interpretation. Furthermore, during preliminary tests we found that introducing the Kullback–Leibler divergence term into the RiboAE loss function added difficulties into the training process.

Thus, the loss function used for the training of the RiboAE, ℓ_AE, is the negative log-likelihood of the input sequence s, where g is the result of applying the RiboEncoder on s and P(s_i | g) is defined in Eq (2).

Evolutionary model

The evolutionary dynamics in our experiment start with an isogenic ancestral population, made of multiple copies of the same SeRANN. However, our framework also supports diverse ancestral populations, with many different SeRANN. We model an asexual haploid population with constant population size N and non-overlapping generations using a Wright-Fisher model. The following steps describe how the next generation of SeRANN individuals is produced from the current generation (Fig 1B).

1. For each individual SeRANN_i, 1 ≤ i ≤ N, the source code s_i is interpreted and trained on a classification task () and a replication task (g → g′). This step models the processes of development and maturation of the phenotype.
2. Each SeRANN_i is given its own genotype g_i and a classification test set of fixed size, {X_j}_j. For each example (X_j,g_i) in the test set, the mature individual SeRANN_i produces a classification prediction and a potential offspring genotype (), where j indexes the test set. Thus, a set of potential offspring genotypes is produced.
   The fertility F_i of SeRANN_i is set to its classification accuracy (the fraction of correctly classified examples out of all classified examples). Its relative fertility is set to . The relative fertility determines the expected contribution of SeRANN to the next generation in step 4, and thus steps 2–4 model the process of natural selection.
3. The number of offspring of SeRANN_i in the next generation N_i is sampled from a multinomial distribution .
4. The genotypes of the next generation are sampled from all the potential offspring genotypes, such that N_i individuals are drawn from the potential offspring of individual i. Steps 4–5 model the process of random genetic drift.
5. The RiboDecoder decodes (translates) each offspring genotype g′ to a new Python source-code string s′. Together, steps 2–6 model the processes of replication and reproduction.

The set of decoded source-code strings s′ is the offspring population, a new population of N SeRANN individuals that replaces the parental population. This offspring population is possibly different from the parental population in terms of genotypes and phenotypes due to replication errors (i.e., mutations), selection and genetic drift. To propagate the next generations, we repeat the training, reproduction, replication, and decoding operations in steps 1–6 (Fig 1).

Experimental setup

We conducted six evolutionary experiments: three with population size N = 100, and three with population size N = 1,000. Our main experiment, with N = 1,000, continued for 6,000 generations. In all experiments, we used genotype length k = 100 and SeRANN source code length shorter or equal to 350 language tokens (see Fig A in S1 Text). Genotype length, population size, and number of generations can all be chosen by the modeler, and were chosen in this study to allow for feasible computation time with the resources available to us: we found that longer genotypes required larger neural networks to replicate them, which required more computation time for training; whereas shorter genotypes were too short to faithfully encode even small neural networks. We executed a short preliminary experiment with multiple different ancestors, selected the ancestor that had the most viable descendants after 100 generations, and used it as the ancestor of the population in the main experiment.

Because of the high computational requirements due to the parallel training of thousands of neural networks, we employed GPUs to execute the simulations. For the experiments with N = 100 we used a single NVIDIA Titan V GPU. For the experiments with N = 1,000, we used a machine with eight NVIDIA Tesla V100 GPUs. An average generation took 9 minutes with N = 100 and 17.5 minutes with N = 1,000. The number of trainable parameters in each SeRANN was limited to 2,000,000 to prevent excessive use of GPU memory. Any SeRANN that exceeded this limit was removed from the population, and its fertility was set to 0. The fraction of such SeRANNs in the population remained very low (about 0.03%), suggesting that the number of parameters did not have a significant effect on networks performance.

For the image classification task that determines the fertility, we used the MNIST handwritten digit image data set [20]. The training set included 57,000 example pairs of handwritten digit image X and a synthetic bit-string genotype g, and SeRANNs were trained for five epochs (each epoch is one complete pass on the entire training set). The evaluation set included 3,000 example pairs that were not used during training, and fertility was determined according to the SeRANN classification accuracy on this evaluation set.

Fitness definition

The fitness of a SeRANN individual is determined by two factors. First, the classification accuracy of SeRANN on an evaluation set determines its absolute fertility, F. Second, some of the SeRANN offspring will not be successfully interpreted as valid Keras neural networks due to replication errors. The offspring survival rate of a SeRANN, V, is thus defined as the probability that its offspring are valid.

Let F_i and V_i be the fertility and offspring survival rate of a SeRANN_i, respectively. Recall that the relative fertility is where j indexes all individuals in the population. The distribution of the number of newborn offspring of SeRANN_i, N_i, is given by a binomial distribution with parameters N (population size) and f_i,

The distribution of the number of surviving offspring of SeRANN_i, , conditioned on the number of offspring, N_i, is given by a binomial distribution with parameters N_i and V_i,

Thus, E[|N_i] = N_iV_i and E[N_i] = N f_i. From the law of total expectation, the expected number of surviving offspring is

Let p_i and p’_i be the frequencies of SeRANN_i in the current and next generation, respectively. Thus, the relative fitness of SeRANN_i is the expected change in frequency of SeRANN_i due to selection, where is the population mean fertility.

However, the relative fitness w_i is unsuitable for comparison of individuals from different generations because of the dependence on the population mean fertility, , which changes over the course of evolution. Therefore, we also use the absolute fitness, (3)

Fertility and offspring survival rate estimation

The training process of a neural network generally starts with random parameter initialization. In SeRANN, this initialization can be interpreted as the stochastic aspect of the development and learning processes, which stems from random environmental differences and chance. Two SeRANNs, completely identical in their architecture, may differ in their parameters (network weights) at the end of training. Thus, to accurately assess SeRANN’s properties, such as its classification accuracy or replication fidelity, a single measurement obtained during the exper- iment is insufficient. Moreover, mutations are rare, the number of offspring of each individual is low (1.59 on average), and the genotype length is short (k = 100), and therefore estimating the mutation rate from offspring produced during the experiment is likely to underestimate the mutation rate.

Therefore, we performed multiple evaluation cycles of SeRANN properties retrospectively. In each evaluation cycle, the individual SeRANN parameters were initialized from a random distribution, and then the individual was trained and evaluated. The initialization and training steps were identical to those performed during the experiment. After training, the individual classification accuracy was evaluated on 3,000 unseen classification examples (as performed during the experiment). Then, the individual SeRANN was given its own genotype paired with random classification examples and the replicated genotypes were examined to estimate the offspring survival rate and the mutation rate.

The offspring survival rate was estimated by first decoding the replicated genotypes to source codes and testing their validity: each source code was executed by the Python interpreter, together with the required supporting code, and was marked valid only if no errors were raised by the interpreter. The offspring survival rate was then estimated as the fraction of offspring of a genotype that were valid. The mutation rate of a genotype g was estimated as the average Hamming distance H(g,g’) between g and its offspring g’ (i.e., the number of sites in which the two genotypes differ), , where N_g is the number offspring of genotype g.

Epistasis estimation

Let g_i be the vector containing the i-th site in the genotypes of h = 126,158 unique genotypes from our main experiment (i = 1,…, k and k = 100). Let F_g be the fertility (classification accuracy) of genotype g. For each pair of sites in genotype g, i, j = 1,…, 100 such that i≠j, we define the linear regression model for predicting F_g as (4) where b, a_i, a_j, and a_i,j are scalar coefficients.

We are especially interested in a_i,j, which estimates the effect of the interaction between sites i and j on fertility. Therefore, we also apply a simpler model in which we fix a_i,j = 0, that is, we assume no interactions between sites i and j. Using the likelihood-ratio test, we test if the evidence supports a non-zero interaction coefficient. Since we perform 9,900 likelihood-ratio tests (one for each pair of sites in the genotype, which has 100 sites), we adjust the p-values using FDR correction [29].

We fitted these models by minimizing the negative log-likelihood with the least-squares method implemented in SciPy. Fig G in S1 Text shows the estimated interaction coefficients for every pair of sites in the SeRANN genotype.

Mutational robustness

The mutational robustness of mutations in genotype g is given by (5) where g′ iterates over all possible mutants of genotype g, M_g→g′ is the number of times genotype g mutated to genotype g′ during the experiment, and W_g is the absolute fitness of the genotype g (Eq 3), based on 50 classification accuracy measurements and 5,000 offspring survival rate evaluations per genotype. Thus, λ(g) is 0 for an optimal genotype (W_g = 1) with only lethal mutants (W_g’ = 0), and 1 for a genotype with only neutral mutations (W_{g’ =} W_g).

Overview

We developed SeRANN, a new neural-network-based computational self-replicator with implicit and endogenous replication and reproductive success and explored its evolutionary dynamics. We observed the emergence of universal evolutionary phenomena: changes in both genetic and phenotypic diversity; adaptation by fixation of both synonymous and non-synonymous beneficial mutations; clonal interference, selective sweeps, and genetic hitchhiking between concurrent mutant alleles; epistatic interactions between genetic loci; and a bi-modal distribution of fitness effects. Furthermore, we observed the evolution of the replication mechanism itself: the evolution of a decreased mutation rate and increased mutational robustness. In the following we summarize these phenomena, with additional details provided in the supporting material.

Adaptation

Adaptive evolution is clearly demonstrated by the increase in the population mean fitness and the decrease in the population mean mutation rate (Fig 2). Indeed, we expect both to improve over time: fitness due to natural selection favoring fitter genotypes; and the mutation rate due to indirect selection favoring genotypes with lower mutation rates [21,30]. Though it decreased, the mutation rate (measured in mutations per genome per generation) remained orders of magnitude higher than that of viruses and microbes [31] and similar to that of microbial populations with mutator alleles [32], suggesting further improvement could be achieved with more time or longer genotypes.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Adaptative evolution.

(A-C) In the first 270 generations, the population mean fitness increased and the mutation rate decreased, whereas the mutational robustness fluctuated. (D-F) Starting around generation 1,350 (dashed line), major adaptations occurred in all three metrics. The number of accumulated mutations increased, but the number of network parameters and the loss_weight value fluctuated rather than trended (Fig H in S1 Text). (G) Population mean fitness (x-axis) vs. mutation rate (y-axis) over time (color, brighter is later). Panels D-G show a rolling average over 10 generations. Fitness: survival rate × fertility (Eq 3); Mutation rate: number of mutations per genotype replication; Robustness: one minus average absolute fitness difference between parent and mutant offspring (Eq 5).

https://doi.org/10.1371/journal.pcbi.1012004.g002

The observed adaptation of the population can be explained by the appearance and fixation of 24 mutant alleles (Fig 3); however, the mutation rate was high enough for the population to remain genetically diverse so that no genotypes fixed during the experiment, and the genetic richness (i.e., number of unique genotypes) remained over 40 (out of 1,000 individuals; Fig E in S1 Text). Of the 24 mutant alleles that fixed, 16 were non-synonymous, that is, they were expressed in architectural and hyper-parameter modifications to the source-code phenotype. The rest were synonymous, with no effect on the source code (Table E in S1 Text), likely the result of genetic hitchhiking [33].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Allele frequency dynamics.

(A) An initially high mutation rate led to early appearance and fixation of 13 mutant alleles in the first 270 generations. (B) Starting from generation 1,350 (dashed line), 11 mutant alleles fixed: four very rapidly, and then roughly one every 585 generations. (C) The time to fixation (y-axis) of 24 mutant alleles that successfully fixed (blue lines in A and B) declines with their fitness relative to the population mean at the time they appeared (x-axis; gray line shows linear regression; Pearson correlation ρ = −0.6, P < 0.002). Allele 38 (pink) is the only one to reach fixation with relative fitness <1, at 0.988. Gray lines denote mutant alleles that never reached 90% frequency. Frequency curves smoothed for visualization using a rolling average over 10 generations. See Fig J and Table E in S1 Text for details on fixed mutant alleles.

https://doi.org/10.1371/journal.pcbi.1012004.g003

As expected, mutant alleles with larger fitness advantage tended to fix faster (Fig 3C). Moreover, longer fixation times occurred during the initial phase of the experiment (first 270 generations), during which genetic richness was high and 13 mutant alleles fixed, including five synonymous mutants, indicating the effects of genetic hitchhiking [33] and strong clonal interference [34].

Additional evidence for clonal interference [34] and for soft selective sweeps [35] is given by the dynamics of the mutant allele at site 62 of the genotype (soft sweeps occur when multiple copies of a beneficial mutation, rather than a single copy, appear and fix together). This mutant allele started spreading in the population at generation 1,486, reaching frequency >90% within 940 generations (black line in Fig 3). Three genotypes carrying this mutant allele reached high frequencies (60%-75%) but none fixed in the population (green, blue, and orange in Fig 4). Although these genotypes could be considered an evolutionary lineage, they appeared more than once during the experiment, from multiple and different parental genotypes, and indeed they were not direct offspring of each other. Rather, they were separated by at least two mutation events, and the intermediate mutant genotypes (i.e., the single mutants; purple, pink, peach, and brown in Fig 4) did not reach high frequencies (<15%), as common during adaptation on rugged fitness landscapes [36], for example due to stochastic tunneling, in which genotypes are part of a successful evolutionary lineage despite never reaching high frequencies [37].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Genotype frequency dynamics: site 62.

Due to clonal interference, the total frequency of all genotypes carrying a site-62 mutant allele (black) is much higher than the frequencies of any site-62 genotypes; three genotypes that reached 60% are in green, blue, and orange. The inset shows the fertility and mutation rate of the different genotypes decreasing over time, although the decrease in fertility (0.14%) is minor compared the decrease in mutation rate (19.43%; Table F in S1 Text). These alleles and genotypes appeared more than once during the experiment, in multiple and different parental genotypes (e.g., m62 occurred 10,702 times during the experiment). Frequency lines smoothed for visualization using a rolling average over 10 generations.

https://doi.org/10.1371/journal.pcbi.1012004.g004

Diversity

We observed 103,666 unique SeRANN phenotypes (source codes). Surprisingly, the number of unique genotypes was approximately 22% higher than the number of phenotypes, indicating that synonymous mutations, which do not affect the phenotype, are common. The number of unique genotypes in each generation consistently decreased during the experiment (Fig E in S1 Text) likely due to (i) a gradual decrease in mutation rates, which caused less mutations to accumulate; and (ii) selection purging unfavorable genotypes from the population.

Evolution of the mutation rate

Overall, we observed a significant reduction in the mutation rate. The start of the experiment (first 15 generations) was characterized by a high mutation rate, between 2 and 7 mutations per genotype replication, on average (Fig 2B), due to a high loss_weight value in the ancestor, which controls the balance between the fertility and replication tasks during SeRANN “training”. Therefore, mutant alleles that reduced the loss_weight value in favor of replication contributed to a decreased mutation rate and can thus be considered anti-mutator alleles, similar to those found in microbial and viral populations [38]: for example, in the long-term evolutionary experiment with Escherichia coli, two MutY mutant alleles that appeared on the background of a MutT mutator allele reduced the mutation rate by ∼56% and ∼36% [32]. Additional SeRANN anti-mutators appeared later, reducing the mutation rate at the cost of decreased fertility (Fig 4 inset, Table F in S1 Text). In two of these cases, the reduced mutation rate can be explained by a 0.002 reduction in the loss_weight value. The resulting decreases in the mutation rate (1.8–13%) were more dramatic than the decreases in fertility (0.01–0.07%), which may explain the evolutionary success of these of anti-mutators. By the end of the experiment (last 500 generations), the population mean mutation rate was as low as 0.04, a 175-fold reduction compared to its initial value. The variance-to-mean ratio of the mutation rate (across the population) also decreased over time, from ~1.5 at the beginning of the experiment to ~1 at its end. The latter is in agreement with the usual assumption that the number of mutations per genome per generation is Poisson distributed.

The mutation rate was not uniform across the genome. A minority of sites received much more mutations, with a mutation rate above 10⁻². Interestingly, the sites with the highest mutation rates were located at the two ends of the linear genome (position 0 and 99), which could be explained by the convolutional layer used in SeRANN to process the genome (Fig M in S1 Text).

Furthermore, we tested the effect of providing random images instead of MNIST images on the evaluated mutation rate. We found that the mutation rate with random images was 2.3-fold higher (paired t-test, P<3×10⁻⁶, Fig N in S1 Text), and that this increase in mutation rate is correlated with classification accuracy (Pearson correlation, ρ = 0.27, P = 0.0087). Therefore, the classification task input (image) affects the performance of the replication task, indicating that the two tasks are dependent, and that unexpected classification inputs can act as “mutagens”.

Surprisingly, we also found a synonymous anti-mutator allele in a genotype that reached relatively high frequency (but did not fix; green line in Fig 4). This allele did not affect the source code, but nevertheless had a significant phenotypic effect, reducing the mutation rate by ~16% (g62-44 in Table F in S1 Text). This synonymous anti-mutator may be compared to the effect of DNA-sequence context on mutation rates [39], such as in CpG hypermutability [40].

We also identified mutator alleles, which increase the mutation rate. For example, a mutant allele that appeared in generation 805 caused a major source-code modification that switched a neural network layer from the replication task to the fertility task (Fig 5). This increased the fertility (from 0.93 to 0.96) but also dramatically increased the mutation rate (from <0.02 to 5.5 mutations per genotype replication). A similar phenomenon was described in Pseudomonas aeruginosa, where mutator alleles (mutS, mutY, and mutM) not only increase the mutation rate but also increase fitness by providing resistance to hydrogen peroxide [41], and in E. coli, where deletion of mutS not only produced a mutator phenotype but also extended to delete part of rpoS, a gene well known for its adaptive potential under stressful conditions [42].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Phenotypic evolution over three generations.

The phenotype of three SeRANN individuals from the same lineage over three consecutive generations: “grandparent”, “parent”, “offspring”. The grandparent had 14 mutant alleles (sites 6, 11, 31, 34, 35, 37, 38, 65, 66, 71, 77, 84, 85, 91) compared to the ancestor; these are the 14 mutant alleles that fixed up to generation 1420 (Table E in S1 Text). The parent was born in generation 805 with an additional mutant allele in site 0. This mutation caused a replication layer (g_layer in blue) to switch to a classification layer (X_layer in red). This required the layer names to change, as well as the layer type to change from Conv1D to Conv2D. This modification in the source code had a strong effect on the fertility (up from 0.93 to 0.96), the mutation rate (from <0.02 to 5.5) and the offspring survival rate (down from 0.9988 to 0.063). Despite the low survival rate, the parent managed to produce an offspring in the next generation (805), with 6 additional mutations. Three were reversions back to the ancestor allele at sites 0, 37, and 65. Three were forward mutations to mutant alleles at sites 22, 47, and 83. These three alleles fixed in the population by generation 1439. The effect of the 6 new mutations was to revert the earlier change that occurred in the parent as well as a 0.002 increase in the loss_weight (higher values favor classification). The offspring fertility, at 0.94, decreased compared to the parent, probably due to the layer switch, but increased compared to the grandparent, probably due to the increase in loss_weight (from 0.1026 to 0.1046). The offspring mutation rate, on the other hand, decreased back to a low level (<0.02) and therefore the offspring survival rate increased to 0.9958, much higher than the parent, though not as high as the grandparent. The offspring genotype was the dominant genotype of the three new mutant alleles (sites 22, 47, and 83, see Table E in S1 Text). See Fig D in S1 Text for source codes.

https://doi.org/10.1371/journal.pcbi.1012004.g005

However, in our case the high mutation rate led to a very low offspring survival rate (0.06, down from 0.99 in the parent). Interestingly, six mutations occurred in one of the mutator offspring, of which three were reverse mutations (i.e., toward the ancestral allele). The combined effect of these six mutations was to revert the major source-code modification described above, as well as to increase the loss_weight value by 0.002 in favor of fertility. The mutation rate thus decreased back to its lower level, although at the cost of an overall 0.02 decrease in fertility (Fig 5).

Distribution of fitness effects and mutational robustness

We accurately measured the fitness effect of every mutation and estimated the distribution of fitness effects of new mutations, or DFE [43]. The estimated distribution is bi-modal and very similar to those estimated in yeast [44] and viruses [45], with high densities of lethal mutations and of neutral or slightly deleterious mutations (Fig 6). We also examined the evolution of the DFE during the experiment. Strikingly, the frequency of lethal mutations decreased over time, whereas the frequency of neutral and slightly deleterious mutations increased (Fig 6). This observed evolution of the DFE is evidence for a phenomenon called “survival of the flattest” [7,46,47]: when the mutation rate is high, as it is in our experiment, the robustness of an organism against lethal mutations is especially important for its survival. Therefore, it is evolutionary advantageous to be in “flat” areas of the fitness landscape, in which mutations have small effects on fitness. Although the fitness in “flat” areas might be sub-optimal, minimizing the effects of mutations on fitness is more advantageous when mutations are common.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Distribution of fitness effects of new mutations.

(A) The distribution of fitness effects (DFE) evolved during the in-silico evolutionary experiment: the frequency of lethal mutations decreased (red line) and the frequency of neutral and slightly deleterious mutations increased (blue line). (B) The DFE in the SeRANN population (blue) is bi-modal with peaks for lethal and near-neutral mutations, similar to a DFE estimated in VSV (red; vesicular stomatitis virus; data from [45]). VSV has more slightly deleterious mutations, whereas SeRANN has more near-lethal mutations, perhaps because the DFE of SeRANN is still evolving. Fitness effects were estimated as the absolute fitness ratio between mutant and parent such that values of 0, <1, 1, and >1 represent lethal, deleterious, neutral, and beneficial mutations. 5,000 parent-mutant pairs were sampled for every window of 1,000 generations for panel A; the average over all generations is shown in panel B.

https://doi.org/10.1371/journal.pcbi.1012004.g006

Further support for “survival of the flattest” is given by the increase in mutational robustness, which is the extent to which a phenotype remains constant in the face of mutations [48,49]. We evaluated the mutational robustness of genotypes that mutated at least once as one minus the expected absolute difference between the parent fitness and the fitness of mutant offspring (see Methods). Overall, the population mean mutational robustness increased over time (Fig 2F) and was positively correlated with the population mean fitness (Fig I in S1 Text). In contrast, Johnson et al. [50] found a negative correlation between fitness and mutational robustness in the yeast Saccharomyces cerevisiae. A possible explanation for our observed positive correlation is that lethal mutations became less likely as the experiment proceeded (Fig 6), which increased mutational robustness, but also increased offspring survival rates and therefore the population mean fitness.

Epistasis

Epistatic interactions between genetic loci were common and had a notable effect on the evolutionary dynamics. This is exemplified by the case of the mutant allele at site 37 of the genotype, m37, which reached fixation in generation 195. Its effect on its dominant genotype (the genotype that carried it and had the largest number of descendants) was a reduction in the loss_weight value in favor of replication. Thus, it can be considered an anti-mutator allele. However, m37 was effectively driven to extinction by generation 1,423 by genotypes that contained the wildtype allele at site 37, w37, and a mutant allele at site 83, m83, becoming the only allele to fix and then become extinct (specifically, the combination m37/m83 appeared in generation 804, and in the next generation back-mutations generated the w37/m83 combination).

We examined all genotypes that carried mutant alleles at sites 37 and/or 83 (Table B in S1 Text). We found that when both mutant alleles are present in the same genotype, the average number of neurons (i.e., units) in the next-to-last layer of the SeRANN decreases by 14–17% compared to the average number of neurons when only one mutant allele is present, respectively. As a result, the average mutation rate increases by 2.5 to 5.5-fold. This has a dramatic effect on offspring survival rate, which drops from 93–98% with one mutant allele to 47% with both mutant alleles. Thus, fitness with both alleles is lower than with either allele, leading to sign epistasis between these alleles [51]. Because the effect of m83 on the mutation rate is stronger than that of m37 and it only requires a minor sacrifice in fertility (Table B in S1 Text), natural selection favors m83, leading to its fixation and to the extinction of m37.

We also performed a global epistasis analysis of all site pairs in the genotype, focusing on fertility, for which we have greater statistical power. We estimated the epistasis coefficient of all pairs of genotype sites, resulting in a distribution similar to those observed in vesicular stomatitis virus (Fig G in S1 Text) and in E. coli [52], with roughly half of the coefficients positive and half negative.

Phenotypic variation

The degree to which the same genotype can develop to different phenotypes under the same environmental conditions is called phenotypic variation. It has been extensively studied in microbial systems, where it can arise from stochastic fluctuations in molecular or cellular processes and attenuated by, e.g., feedback regulation [53]. In SeRANNs, network parameters are randomly initialized before training and therefore multiple SeRANNs with the same architecture trained on the exact same sequence of examples may have different parameters and therefore may differ in their fertility and replication performance.

We measured the variation in SeRANN fertility and offspring survival using repeated cycles of initialization-training-evaluation. We found a wide range of fertility variation (between 0.001 and 0.02; 99^th percentile <0.012) and an even wider range of offspring survival rate variation (between 0 and 0.4; median 0.09). Fertility variation positively correlated with average absolute fitness (Pearson correlation, ρ = 0.52; permutation test, P<10⁻⁴; Fig F in S1 Text). In contrast, the relationship between variation in offspring survival rate and absolute fitness is more complex. Fitness cannot be close to the neither 0 nor 1 when offspring survival rate variation is high, because fitness is estimated by the average offspring survival rate times the average fertility across all evaluation cycles (Eq 3). Thus, despite a lack of correlation between the standard deviation of offspring survival and absolute fitness (ρ = 0.002, P>0.73), the fittest SeRANN individuals have low variation in offspring survival (Fig F in S1 Text).