15—
Some Elementary Attempts at Numerical Modeling of Problems Concerning Rates of Evolutionary Processes:
With R. Schrandt (LA-4573-MS, December 1970)

This report is an attempt to study mathematically models of evolution showing the qualitative difference between the rates of development in sexual and non-sexual processes. A further number of papers were stimulated by this report.

See also the earlier Proceedings of a meeting at the Wistar Institute on "Mathematical Challenges to the Neo-Darwinian Interpretation of Evolution" held in April 1966 where a paper of mine called "How to formulate mathematically problems of rate of evolution" was published. (Wistar Institute's Symposium Monograph No. 5, pp. 21-33, 1967.)

The work of the above report, done with Schrandt, was performed before the Wistar meeting. (Author's note). *

Abstract

An account of numerical work prepared on electronic computers-the problems concerned the ratio at which favorable mutations spread throughout the population

* This report is also reprinted in the proceedings of a Los Alamos conference dedicated to Stan Ulam, Evolution, Games and Learning that occurred May 20-24, 1985 in Los Alamos. (Eds.)

subject to a "survival of the fittest" mechanism. Models of asexual reproduction showed the expected linear growth in the number of improvements. The bisexual process greatly accelerated the average acquisition rate.

I—
Introduction

In this report, we shall present an abbreviated account of calculations performed by us in the mid 1960's. These calculations were preliminary and intended merely as the zeroth approximation to the problem concerning rates of evolution-a process which we have here severely stylized and enormously oversimplified. A mention of the results of such calculations in progress at that time was made at a meeting in 1966 at the Wistar Institute in Philadelphia by one of us. The discussion there, as reported in the proceedings of the meeting, was rather frequently misunderstood and the impression might have been left that the results somehow make it extremely improbable that the standard version of the survival-of-the-fittest mechanism leads to much too slow a progress. What was really intended was indications from our computations-simple minded as they were-that a process involving only mitosis, in absence of sexual reproduction, would be indeed much too slow. However, and most biologists realize it anyway, the Darwinian mechanism together with mixing of genes accelerate enormously the rate of acquiring new "favorable" characteristics and leave the possibility of sufficiency of the orthodox ideas quite open. Numerous requests addressed to us for the elucidations and details of the numerical setup made us decide to give this account of our computations.

Perhaps the greatest uncertainties-the strongest objections to any calculation of the sort described in the pages that follow-must concern the values of the constants which are assumed initially or should indeed concern even their meaning in the interpretation we have chosen. We have tried to interpret the survivability of individuals by changes in the number of offspring which carry the species in time measured by a discrete succession of generations. The value of "favorable" mutations was mirrored in the increased proportion of offspring. The same, needless to say, is true of the frequency of favorable mutations. We have disregarded the lethal and the unfavorable mutations. We assumed a special form of the advantages which an individual holds relative to the rest of the population by comparing the number of his "improvements" with the average number present at that time in the

population. We assumed a proportionality law, again arbitrarily. In some problems we have penalized the individuals whose score in the improvements was less than the average; in some of the problems we considered only the positive excess as leading to a greater number of offspring. Another debatable procedure is the way we have handled the growth of the population by normalizing periodically the total number of individuals to a constant figure. If the number of individuals holding a certain number of "advantages" after normalization dropped below 1, we summarily dismissed such representatives. It should be stressed here strongly that this procedure makes it very hard to find an analytic model equivalent to our numerical work. We do not have any clear idea of the necessary scaling laws concerning the effect of changing the constants alpha, kappa (defined below) and the size of the population. All this is true throughout all the problems. In the calculations involving combinations of genes from both parents, further assumptions were made of independent inheritance of the "improved" genes, etc. As will be seen in the description of the individual problems, we have chosen successively less unrealistic assumptions. Clearly in the counting of the new improved genes coming from the "father" and from the "mother", one has to take care not to count the same "improvements" from each twice. As will be seen, this precaution has the effect of slowing down the at first seemingly exponential growth

into something more like a quadratic function of time (we have studied throughout the calculations the number of "improvements" present in the population as a function of time, that is to say of the generation number).

In order to get a feeling for the dependence of the results on the values of the constants, more such computations must be tried in the future and additional variables have to be considered-certainly the "kind" of the improved or favorable new gene has to be taken into account. A most important question concerns the existence of new genetic instructions involving perhaps logical prescriptions, that is to say recipes for operations and actions of the components rather than merely their chemical composition. An improvement in programming or interpretation of action by a gene or group of genes may be equivalent to a very large number of the "favorable changes" with which our computations have dealt so far.

Our first problem with the code name ADAM concerned asexual reproduction. We feel that the time scale to acquire a characteristic in an organism, such as the development of an eye, by a sequence of consecutive favorable mutations, is extremely long if one does not resort to something like sexual mating in the population. In the following rough and elementary estimate, the constants assumed are crude, but err toward "faster" evolution than what is to be expected.

Definitions: Let T = time of existing life (~10⁹ yrs.); T = time for one generation say (-3 days); G = the number of generations=T/T=1011; N= the existing population size (-10¹¹ ); K = the total number of "favorable mutations" necessary to produce the desired characteristic (--106); a = the chance of a favorable mutation per individual per generation (~10-1°); y = the "value" of a single favorable mutation expressed as a survival rate (~10-6). That is, an individual having this mutation would have (in expectation) (k + 7y) descendants, versus k descendants for an individual not having this mutation.

Therefore, in the first generation, the expected value of the population that could have one mutation is N ca = 10. In 1/7y =106 generations, a sizeable portion (approximately l/e) of the population would have this mutation, and in about 10⁷ generations, most of the population would have it. But the time to acquire all the mutations would be about K. 107, or 1013 generations, which is like the age of the universe.

II—
ADAM

In order to study, on a computer, the rate at which a population can acquire a sequence of mutations we needed a set of more amenable parameters, which, it is hoped, could eventually be scaled down to the set given above. For the first problem, called ADAM, we used the following set: N = 100, K = 100, a = .02, and 7 = .1.

The method was as follows: In any one generation, each member of the population N had a probability a of acquiring one new independent mutation. Each individual then had one child, with a probability of extra children given by y(Kg - Ko), where Kg = the total number of different mutations possessed by this individual and, Ko = the minimum number of mutations possessed by any individual in the population. (If n/-y < Kg - K <(n + l)/7y, n = 1, 2, 3, ..., the individual had n extra children, and a probability of y(Kg - Ko- n/y) for an (n+ l)^th child.)

The children were then given the number of mutations possessed by their parent. The parent population then was assumed to have died, and the children formed the new population. The numbers of

mutations in the population were recorded in categories by counts ni, i = 0, 1, 2, 3, ..., where ni = the number of individuals having a total of i mutations, and Eni =N, the size of the new population. In the next generation, each member of this new population had the same probability a of acquiring another new mutation, and had children according to the above recipe. These children with their number of mutations recorded then became the population for the succeeding generation, etc.

It was necessary to renorm the population periodically, since the number of children increased in each generation. This was done by reducing the count ni in each category by 1/2, to the nearest smaller integer, when the population reached 200, which is double the initial population.

The weighted average number p of mutations possessed by the population was computed for each generation from the categories ni. This average was then plotted as ordinate against the generation time as abscissa. The slope of this curve is then the rate at which the population can acquire a sequence of mutations, as a function of the parameters a and y.

This rate of acquiring mutations turned out indeed to be linear. For the parameters a = .02 and y = .1, the slope was about .1; or a majority of the population acquired an additional mutation every 10 generations. Several problems were run with smaller values of y, that is, y' = f -y, where f = a fraction. The graphs were all linear, with decreasing slopes, which decreased more closely with Vf than with f itself. There was no appreciable change in the slope by doubling the initial population to 200. Figure 1 shows a plot of 3 cases: N = 100, a = .01, and y = .1, .05 and .01 respectively.

A second version of ADAM was run with the Ko in the probability recipe defined as the average number of mutations in the population, instead of the minimum number. Those individuals having fewer than the average number of mutations (Kg < Ko) had a probability of 7(Ko - Kg) of no children, (if Ko- Kg> 1/7, they had no children deterministically), and a probability of 1 - 7(Ko - Kg) of one child. The individuals with more than the average number of mutations had the same probability for extra children as the one defined above. This version required fewer renormings of the population and it led to a somewhat greater slope than the Ko = min. recipe. The graph was still linear. Figure 2 illustrates the two versions with the same parameters N = 100, a = .02, and - =.1. The recipe with Ko = average was used in all subsequent problems.

NUMBER OF GENERATIONS Fig. 1.

III—
EVE

In the second class of problems, we introduced reproduction in the population, and also, what seems important, fluctuations in the manner in which the offspring received mutations from their parents. This problem was naturally called EVE. The initial parameters used were the same as in ADAM. The population acquired new mutations according to the probability a, as before. A random mating of the population N was then defined, resulting in N/2 pairs of individuals. (For N odd, the population was arbitrarily reduced by one to obtain (N - 1)/2 pairs.) Each pair then constituted a set of parents. The number of children from each set was again determined by the probability function a(Kg - Ko). Here Kg = the total number of mutations possessed by the set of parents, and Ko = the average total number of mutations of all the pairs of this mating.

* K, MINIMUM ^ K. AVERAGE NUMBER OF GENERATIONS Fig. 2.

The children were produced in pairs. Thus if a set of parents had Kg > K₀ , they had two children for certain, and with a probability given by y(Kg - Ko - n/y), n= 1, 2, 3, ..., where n/y < Kg - K < (n + l)/a), of 2n extra children. If for the set of parents, Kg < Ko, they had, with the probability y(Ko - Kg) no children, and with the probability 1 - y(Ko - Kg) two. Again if K - Kg> l1/, they had no children.

The number of mutations acquired by each child was obtained under a binomial distribution centered about the average number of mutations of the parents. If the parents had a count of (x + y) mutations, the number for each child was obtained under the distribution centered at (x + y)/2, with its minimum at zero and maximum at (x + y). The number of mutations for each child separately was determined under this distribution. Thus a child could possibly obtain as much as the sum of numbers of mutations of its parents. This recipe

involving fluctuations in the inheritance would, we thought, speed up the rate of acquisition of mutations (compared to always giving the offspring (x+y)/2 favorable mutations).

The parents having died, the children became the new population. As before, they were classified in terms of counts of the number of individuals having each number of mutations. As before, the individuals had the chance a new mutations, and were mated at random in N/2 pairs. These pairs then had children whose mutations were again determined from the binomial distribution, and these children constituted the population for the following generation, etc.

In our random mating, the sex of the individual was not distinguishable. No attempt was made to keep members of the same "family" from mating. Their number of mutations was not necessarily the same, since it was determined separately for each child under the probability distribution. (The identity of the family was lost once the members were classified according to their mutations.)

The norming of the increasing population was done in the same way as in ADAM: All categories were halved when the population doubled.

The rate of acquiring mutations turned out to be much faster than in ADAM, and appeared to be exponential. In Fig. 3 we have plotted this rate on a semilog scale for four problems: N = 100, a = .02, and y = .1, .05, .025, and .01 respectively. The reduction in acquisition rate with y was somewhat similar to that in ADAM.

In these problems no attempt was made to keep the histories of the different mutations. We define a to be the rate of acquiring new mutations, but we divide the population only in terms of numbers of individuals having a fixed number of mutations. The children too acquire mutations only under the distribution of the total count of the parent's mutations. Thus one might suspect that the exponential rate of acquisition could be due to a doubling of the identical mutations possessed by both parents.

A—
EVE-PQ

To correct for this, we computed an expected number v of mutations that the parents should have in common. This is given by v = (T₁ T₂)/S where T1 = the total number of mutations possessed by one parent, T2 = the total number of mutations possessed by the other parent, and S = the total number acquired by the entire population. (The

NUMBER OF GENERATIONS Fig. 3.

total S is accumulated as each individual in each generation has the probability a of acquiring a new mutation. If Ns is the average population in each generation, after k generations, Sis approximately N . k · a.)

We allotted this number v to the children for certain and then played a game of chance for additional mutations using the reduced binomial distribution centered at the midpoint of the total count of independent mutations possessed by the parents. In this manner, we count more correctly, the number in common only once.

This method leads to a slower rate of acquisition; it is still exponential in the beginning but tails off to something like a quadratic function, as more mutations must of necessity be held in common. Results from four problems are plotted in Fig. 4 on a semilog scale. The

NUMBER OF GENERATIONS Fig. 4.

cases PQ1 and PQ2 have N = 100, a = .02 and y = .1 and .05 respectively. The case PQ3 is the same as PQ2 except that the sample size was doubled, N = 200. It had a somewhat higher value for p,, but the slopes are like those of PQ2. The last case, PQ4, had the same parameters as PQ3 except that a was cut to .01. This problem was run to 290 generations and showed a definite bending over of the curve to almost a linear rate of acquisition after 230 generations as the population had more and more mutations in common.

B—
EVE-PM

In the two versions of EVE already discussed, the population mated at random uniformly. We next considered a version where the mating was not uniform, i.e., preferential in the following sense: We arbitrarily

divided the population equally into three groups ranked according to their number of mutations, i.e., the first had the individuals with the greatest number of mutations. We specified that 3/4 of the population in each group mate at random uniformly within their own group. The remaining 1/4 would mate at random with equal probability from either of the two remaining groups. For example, if we name the groups A, B, and C, an individual from group A would have a 3/4 chance of acquiring a mate from group A, and a 1/8 chance of a mate from each of the groups B and C. We called this problem EVE-PM. The EVE-PQ method was used to estimate the mutations in common and to count them only once.

The rate of acquisition of mutations should reflect this preferential mating. A comparison of the curve PM1 of Fig. 5 with PQ1 of Fig. 4, (both with the same parameters), shows indeed that the initial acquisition rate of mutations under preferential mating is much higher than under random mating. But the curve PM1 tails off very rapidly after 100 generations, and at 200 generations the mutation rate is almost the same for the two problems. This indicates that at this point most individuals have acquired most of the mutations available in the total population, so that preferential mating has the same effect as the uniform mating. Our small sample size of 100 in part causes this phenomenon to occur so soon.

The computing time goes up rapidly with the initial population; one problem of preferential mating with an initial sample of 400 was run to 96 generations. The result is shown in the curve PM2 of Fig. 5, with the same parameters as PM1 (except for the population size). The preferential mating has a greater initial effect in the larger population, but the slope of this curve too is beginning to decrease.

The mutation acquisition rate in all the EVE problems of sexual reproduction can apparently be divided into three stages: * An initial exponential rate, as few mutations are held in common (compared to the size of the population). * A rate, roughly quadratic, as more mutations are held in common by the parents. This number in common is approximated by computing the expected intersection of the number possessed by each parent, assuming the parents had acquired them independently. * A terminal rate, almost linear, as most of the mutations in the population are in common. If all the mutations were in common, the subsequent rate of acquisition must be linear, since new mutations are obtained only by the a rate of new acquisitions, which is a linear function.

NUMBER OF GENERATIONS Fig. 5.

C—
EVE-POS

In order to check these assumptions somewhat, a problem was run in which histories were kept of each new mutation. This was done by representing each new mutation as a bit position in a matrix of words in the computer. Each individual in the population and its children had their mutations recorded in such a matrix, that is to say mutation was specifically identified.

We called this problem EVE-POS. The sexual reproduction scheme was the same as before. The population was mated at random, and each set of parents had a probability of having extra children given by y(Kg - Ko). In both Kg and Ko the actual mutations in common were counted only once. This number was known for each mating, so no approximation for the expected intersection was needed.

If both parents possessed the i-th mutation, it was then given to each of the children. If neither parent possessed this mutation, it was not acquired by any of their offspring. If only one parent possessed this i-th mutation, each child had a probability of 1/2 of acquiring it.

The norming for this problem was different from before, namely, when the population doubled, each individual was given a half chance of surviving.

With this recipe for receiving mutations, any individual mutation can be lost to the population, since the parents die off in each generation. For example, the k-th mutation is initially acquired in the a recipe by one individual. If this individual and its mate have no children, the mutation is lost. If they have n children, there is a probability of (1/2)" that none of the children get it, in which case it is lost in the next generation. There is a chance that this particular mutation will be lost in each subsequent generation, although these probabilities are getting smaller. The k-th mutation is initially acquired only once, and by only one individual. The mutations that survived were "packed" in the bit positions of the matrices. This relieved the space limitation in the memory and allowed the problem to be run much further in time.

It was discovered that approximately 80% of the total of new mutations acquired by individuals were "lost" after matings in subsequent generations. Thus the parameter a in problem EVE-POS has a different meaning from that in EVE-PQ. In EVE-POS it denotes the probability of an individual acquiring a new mutation. In EVE-PQ it denotes the probability of acquiring a mutation that survives and will eventually be acquired by the entire population.

Figure 6 shows the mutation rate for two cases of EVE-POS, plotted on a log scale. The parameters are a = .02, y = .1, and N = 100 and 400, respectively. The case POS1 (N = 100) was run to 251 generations, and POS2 (N = 400) to 90 generations.

Note that the acquisition rate and the values of p are considerably smaller in POS1 than in PQ1, which has the same parameters. This is because of the different interpretation of the parameter a. A reasonable comparison with POS1 would be to run PQ1 and compute our v= (Ti • T ₂)/S', where S' = .2S, since about 80% of the mutations are lost.

The acquisition rate for POS1 becomes linear after about 100 generations. The relatively small sample size is a contributing factor. Some statistics were compiled on the distribution of the available mutations. They are given in Table I.

An estimate was made of the expected number of mutations held in common compared to the actual number held by an average set of parents. The expected number was computed assuming that each

NUMBER OF GENERATIONS Fig. 6.

parent had the average number p of mutations. Then v =p² /S', where S' = the total number of surviving mutations. The actual number held in common was about 1.3 times this expected number, after 251 generations.

In the problem with the larger sample size of 400 (POS2), the mutation acquisition rate remains approximately quadratic through 90 generations. At this point 18% of the mutations were held by at least 50% of the population, versus 36% for the population of 100. For this larger sample size, at 90 generations, separate distributions were kept of mutations acquired in the first 45 generations, versus those acquired in the last 45 generations. This data and some statistics on those mutations held in common are given in Table II.

15Elementary Attempts at Numerical Modeling of Problems of Evolution

TABLE I Problem-EVE POS1-N =100, a = .02,y = .1 After 90 generations, there were 55 surviving mutations out of 266 acquired mutations: 79% were lost. Distribution of the 55 mutations: Min. - 15 (least number held by any individual) Aver.- 21 (average number held-this is the number p that is plotted) Max. - 27 (greatest number held by an individual) 36% of the mutations were held by at least 50% of the population 1.8% of the mutations were held by the entire population.

After 251 generations, there were 130 surviving mutations out of 744 acquired: 85% were lost. Distribution of the 130 mutations: Min. - 87 Aver. - 94 Max. -101 74% of the mutations were held by at least 50% of the population 41% of the mutations were held by the entire population.

The actual number of mutations held in common was approximately 1.3 times the expected number.

TABLE II Problem-EVEPOS - N400,a .02, y. After generations, there were surviving mutations out of acquired: were lost. Distribution of the total mutations: Min. - Aver. - Max. - Distribution of mutations initially acquired in the first generations (: Min. - Aver. - Max. - Distribution of mutations acquired in the last generations (: Min. -Aver. - Max. -There were of the mutations held by of the population, held by of the population, and held by at least of the population. At generations, the actual number of mutations held in common was ap- proximately times the expected number.

The above problems give some indication of the rate at which a population can acquire mutations in terms of a finite and rather small sample size, and in terms of the procedures which we adopted for acquiring and transmitting of the mutations. These methods involved using rather large values of the parameters alpha and gamma. The problem remains to find scaling factors. The computing time on the old IBM-7094 for the problem POS2 was over one hour, for 90 generations of growth. The computing time for the PQ code was much faster, in the order of minutes, but the parameter a was in effect much larger than the one in the POS code.

IV—
Summary

The problem ADAM with asexual reproduction gave a linear rate of acquisition of mutations. In reducing the parameter y to f y, where f is a fraction, the acquisition rate was reduced by a factor more like f/7 than like f itself.

In problem EVE with sexual reproduction, the acquisition rate appeared to be exponential if the initial population were large enough. But with a small population, more of the same mutations were held in common by the parents. This caused the rate to change from an exponential to a "quadratic" and eventually to a linear one when most mutations were common to the majority of the population. The problem EVE-PQ involved approximating this number (our formula for v).

The advantage of preferential mating over random mating gave an initially pronounced increase in the acquisition rate, but this was soon offset by the smallness of the population. In effect, as more mutations were held in common, the range of the distribution of mutations became narrow. After that the preferential mating was not much different from the uniform one.

The EVE-POS problem (where we kept a history of the mutations) gave us a measure of the distribution of mutations as a function of their age. It showed that most of the mutations initially acquired by one individual were lost in subsequent matings. This caused a redefinition of the probability a in computing the expected number of mutations held in common. This problem also showed that the actual number held in common was greater than the estimated number v, by a factor of about 1.3 for the sample size of 100, and 2.4 for the sample size of 400. This is not too surprising, since the expected size of the intersection assumes independent sampling, whereas the mutations are acquired by something more like a Markov process.