Simulated self-organisation of death by inherited mutations

An agent-based computer simulation of death by inheritable mutations in a changing environment shows a maximal population, or avoids extinction, at so intermediate mutation rate of the individuals. Thus death seems needed to al for evolution of the fittest, as required by a changing environment.


Introduction
More than a century ago Weissmann argued that ageing and death are needed to make place for our children; and children are in turn needed to allow for Darwinian evolution through survival of the fittest. Kirkwood [1] summarized this and many other theories of ageing, and specific computer models of ageing and death are reviewed e.g. in [2,3], for example the Penna model (as also reviewed in [4]). (A mathematical argument against immortality was given in this sense in [5].) Now we want to understand the need for death through Monte Carlo simulations of individuals. We distinguish between newborns and adults, and take into account environmental changes. They may come from climate change, like ice ages and warmer periods during the existence of homo sapiens. Or they may be caused by migrations of people from one environment to another. A single such environmental change was already used to justify sexual over asexual reproduction [6]. Thus we vary the mutation rate of individuals to find its optimal value. Here "optimum" either means a maximum of the population in a fixed environmental carrying capacity, or survival instead of extinction, depending on which of our two models A and B we use.
In our two models A and B, using sexual reproduction, the genome is represented by two strings of L bits each. They represent the L most serious genetic diseases. Each mutation damaging the phenotype (i.e the health of the individual) reduces the survival probability per iteration by a factor x.
As genetic load we count those bit positions which differ from an ideal bitstring. The latter is initially zeroed, but changes at each iteration with a probability p at one randomly selected bit position, and thus represents the changing environment. For reproduction the two bit-strings of the father are crossed-over at one randomly selected bit position, the same happens for the mother, and then one of the two resulting bit-strings from the father (the gamete) is combined with one of the two from the mother to give the child genome. Mutations are thus inherited from the parents, and m new mutations are introduced at birth to each gamete (if m ≥ 1; for m < 1, one new mutation is added with probability m). N is called the genetic load; more precisely it is the number of loci (bit positions) where the genome is not adapted to the current environment. All changes in the individuals and the environmental bit-strings are reversible.
Model A, discussed first, uses a varying population and finds as an optimal m that mutation rate for which the equilibrium population reaches a maximum. Model B, discussed thereafter, follows a tradition of theoretical biology and keeps the population constant except if during one iteration all adults die out; then we check which mutation rate avoids extinction of the whole population. Further details of the two models will be discussed in the corresponding sections.

Model A with changing population
In model A, each of the individuals survives the next time step (iteration involving all survivors) with probability x N (1 − P/K) where P is the current total population and K is a fixed input parameter, sometimes called the carrying capacity representing limitations (due to lack of food and space) for the growth of the population. Here N bits are not adapted to the current environment. (The Verhulst factor 1 − P/K applies to all individuals, differently from [7].) Our mutations are recessive, which is defined differently in two different versions A1 and A2 of model A. For A1 we take the logical and of the two  bit-strings of the individual, and then count as N the number of bit positions where this logical and differs from the current ideal bit-string. This version is close to [8], and means that heterozygous loci do not count for the genetic load N before the first environmental change, but count after it. For A2 we count for N only those positions where both individual bit-strings agree with each other (homozygous loci) and disagree with the ideal bitstring.
Our K is mostly 2 million, the initial population is K/5, and the resulting equilibrium population is mostly of the order of one million if it does not die out. The two individual bit-strings are mutated independently, each with mutation rate m. Each surviving adult at each iteration gives birth to B babies, which become adult at the next iteration; we used B = 4. Mostly 10,000 iterations were made (100,000 for most cases with m < 0.001), and averages were taken from the second half of this time interval.
Case A1: Figure 1 shows our main result: The population P has a maximum as a function of m at some intermediate m value. Thus neither very small m ("eugenics") nor very large m ("instability") are optimal; an intermediate mutation rate leads to the largest P or the lowest < N >, but also in reality to a finite lifespan.
Instead of applying the Verhulst deaths to all ages, Fig.2 shows the correlation between genetic deaths only and genetic load, by applying the Verhulst deaths only to the births [7]. Data are taken from averages calculated at different time steps of a single run for each value of x.
Case A2: The modified recessiveness defined above for model A2 reduces  The model treats heterozygous loci in the same way as model A2, that is, they never contribute to the genetic load. A slightly different version of this model, in which x was recalculated at each time step to keep constant the fraction of deaths, was presented in Ref. [9].
The results for this model shown in Fig. 4 should be compared to our Fig.  1. We keep the mutation rate of the environment p = 0.01 and the selection strength x = 0.98 fixed and compute the average genetic load < N >, the fraction of the population that dies per time step, and the fraction of perfect, or ideal, genomes in the population as the mutation rate m is varied. We find that there is an intermediate range of values of the mutation rate for which both the genetic load and the death rate go through minima, while the fraction of perfects reaches a maximum. This result matches what was found, in similar situations, in our model A1.
Our main result refers to the need for a strong selection mechanism as a means to enforce a small genetic load: death of the least adapted individuals makes way to fitter ones. In Fig. 5 we show the time evolution of the average genetic load of the population for four different sets of parameters. In all four, we simulate a population of 1000 individuals, each represented by two bitstrings of size 3200 bits each, with a mutation rate at birth of m = 1.0. In case (a), x = 0.98 (weak selection) and p = 0, the environment does not change. The average genetic load starts at 0 (ideal individuals) and grows to a small value of order 1. The distributions of genetic loads are shown in Fig. 6, averaged after the initial 5000 time steps. For case (a), it has a peak at 0, meaning that a majority of the population carries no genetic load, with a small width. When the environment changes with probability p = 0.01 at each time step (case (b)), the average genetic load increases to a value of order 10 and its distribution peaks at a small non-zero value of the same order. Further increase in the rate of environment change to p = 0.02 leads the population to extinction (case (c)). The average genetic load increases rapidly and its distribution widens: it is shown in Fig. 6. The genetic load accumulates thanks to the joint effects of the mutation rate at birth and a fast environment change and, even with a weak selection, leads eventually to extinction. The need for a strong selection is now shown: for the same parameters (p and m) but smaller x = 0.95 (case (d)), the population resists and the distribution of genetic load is very similar to the one in case (b).
The same qualitative results were obtained for haploid asexual populations, but extinction is avoided only for larger populations and stronger selection pressure, similar to [10]. Extinction can then be correlated to features of the distribution of genetic load. It is avoided as long as the average genetic load is not much larger than the width of the distribution. This is more clearly shown in Fig. 7 is the outcome of the simulation. In the same plot we also show the fraction of individuals that die (for genetic reasons only in this model) at each time step. As p is increased, survival of the population becomes more difficult and causes this fraction to be ever increasing.

Conclusion
In all our models, the genetic heritage of a diploid individual is represented by a pair of bit-strings, which undergo mutations at birth, while the ideal phenotype is mapped into a single bit-string. Environmental change is translated into a mutation of this ideal phenotype. The genetic load of an individual is determined by a comparison between its genetic strains and the ideal phenotype. This genetic load determines the death probability of each individual.
Our results come from simulations with a fixed rate of environment change and a fixed value for the parameter that measures selection strength x. We show that population fitness, determined by its size, reaches a broad maximum, while the average genetic load reaches a minimum, for some intermediate range of the mutation rate at birth (model A1). So, nature has self-organised its cellular error correction machinery to ensure a mutation rate within some range.
On the other hand, when the rate of environment change increases, our results are consistent with the interpretation that selection has to get stronger to avoid population extinction (model B).
A more realistic approach would be perhaps to assign a different selective value for each different bit position, since different inherited diseases differ in their danger to survival. However, that modification would introduce so many free parameters that the model would lose its value.    Values correspond to averages taken after 5000 initial time steps, up to 10 6 . Similar results are obtained for L ≤ 32768.