Physica A 252 (1998) 325–335
A model of population dynamics Andrzej Pekalski ∗ Institute of Theoretical Physics, University of Wroclaw, Pl. Maxa Borna 9, 50-204 Wroclaw, Poland Received 24 October 1997
Abstract A model describing the evolution of a sexual population in a selective environment is presented. The population is composed of individuals each characterized by its phenotype and age. Within the standard Monte Carlo simulation technique, we calculate the time dependence of the average tness, average adaptation to the environment and the distribution of the phenotypes. We show that the former quantities exhibit damped oscillations, meaning that in the absence of mutations the population is driven by the selection to a homogeneous one. The distribution of the phenotypes, as well as the adaptation is Gaussian at each time step of the process. The role of the maximum age and length of the genetic strings on the dynamics of the population is also c 1998 Elsevier Science B.V. All rights reserved. discussed. PACS: 0.2.70.Lq; 87.10.+e Keywords: Biological evolution; Monte Carlo simulations; Population dynamics
1. Introduction The problem of describing dynamics of a biological population has attracted attention of physicists in recent years [1,2]. In some cases [3], sexual reproduction takes place at a genotype level and from it a phenotype for an individual is constructed. In a recent paper, Doebeli [4] presented interesting models of quantitative genetics and population dynamics. Out of several problems he studied, we select just one – the eect of varying character on the dynamics of population with sexual reproduction. The basic assumptions of the model are the following. A population of haploidal individuals is considered, where each is characterized completely by the value (h) of its phenotype (i.e. character). Genetic strings have c loci, each of which can be in two states, called alleles, and denoted here as 0 and 1. In other words, it is a string of c bits. The value of the character for an individual is then simply the number of 1’s in its genetic ∗
E-mail:
[email protected].
c 1998 Elsevier Science B.V. All rights reserved 0378-4371/98/$19.00 Copyright PII S 0 3 7 8 - 4 3 7 1 ( 9 7 ) 0 0 6 3 9 - 0
326
A. Pekalski / Physica A 252 (1998) 325–335
string. Genetic reshuing is introduced as a rule according to which an ospring takes an allele (bit value) at each locus with equal probability from both parents. Dynamic rule determining the phenotypic frequency distribution in the next generation depends on the distribution of overlaps of the parents’ phenotypes (assumed binomial) and the tness function. The latter is taken in a simpler [5] and then in a more complicated [6] form. Both are, however, phenomenological. A numerically solved dierence equation shows that the average tness, de ned as the ratio of population size in two subsequent time steps, has an oscillatory character. The period of oscillations decreases with the number of loci. Doebeli also contests common assumption made by biologists (see e.g. Ref. [6]) that the distribution of a continuous character is normal in each generation. Indeed, his results show, for sexual reproduction, a phenotype distribution uctuating between two maxima. Since the analysis in Ref. [4] is carried out on a mean eld level, as he is dealing only with distributions of the character, it seems interesting to construct a similar model but operating on individuals. The problems which seem worthy to investigate are – how robust is the oscillatory character of the average tness time dependence on the details of the model; is the double-peak distribution of phenotypes an artifact; and what will happen if the population is put into an environment which introduces selection acting at the individuals’ level.
2. The model Our investigation will be based on the standard Monte Carlo simulation technique. Consider a square lattice L × L, where initially N (t = 0) haploidal individuals are located. Each of them is characterized by its location (j) at a lattice site, its genetic string (fj ), of c loci, each with two alleles manifesting themselves as either a 0 or a 1, and the individual’s age (wj ). The phenotype of an individual is given by an integer number between 0 and c equal to the sum of the values of the alleles at the c loci. The population is in an environment favoring the 1’s. Hence, a 0 in an individual’s genetic string corresponds to a bad gene and 1 to a good gene. The initial distribution of the phenotypes is random. The algorithm (dynamical rule) is given by the following steps. (1) An individual is randomly chosen. (2) Its adaptation, i.e., the fraction of 1’s in its genetic string of fj is calculated aj =
c 1X fj : c
(1)
=1
Then a random number r ∈ [0; 1] is generated. If r ¿ aj the individual is removed from the system. If not (3) the individual is moved to an empty nearest-neighboring site. If none is available, return to 1.
A. Pekalski / Physica A 252 (1998) 325–335
327
(4) After the move, a partner from one of the nearest-neighboring sites is chosen randomly. If all sites are empty, return to 1. (5) Survival of the partner is checked (Eq. (1)). If successful (6) q osprings are created. Independently for every one of them at each locus an allele is chosen, with equal probability, from the alleles of the parents (genetic shuing). The osprings are located on empty sites inside a square LG × LG centered at the rst parent location. Since in our model individuals may mate after reaching age equal one, putting the osprings close to the parents increases chances of their future mating. It is therefore, in some sense, equivalent to inbreeding in biological systems. If less than q sites are available, there are less osprings. (7) After completing the procedure as many times as there are actually individuals in the population (one Monte Carlo step – MCS), the age of all members of the population is augmented by one. Individuals whose age exceeds the assumed maximum age (M) are removed from the system. The procedure is repeated with dierent random distributions of the phenotypes of the initial population and then the results are averaged. This is a simpli ed version of the model presented already elsewhere [7], where a colonization of an empty ecological niche by a diploidal sexual population has been considered. The parameters of the model are the following: length of the genetic string (c), maximum age (M ), maximum number (q) of osprings, linear size (LG ) of the area inside which the osprings are located and the density of the initial population x(0) = N (0)= (L × L). Several of these parameters have direct relation to the parameters used in Ref. [4], like the number of loci, number of osprings (q) which corresponds to intrinsic population growth . pj is equivalent to individuals’ ability to cope with the environment. Since in Ref. [4] one deals with characters distributions, osprings are located over the whole system, hence there LG = L. An important factor in Ref. [4] is the parameter re ecting competitive interactions among dierent types of the phenotypes. In our model the competition is brought at the level of individuals – with the limited carrying capacity of the system ( nite size of the lattice with hard boundary conditions), the reproduction rate depends on the adaptation of parents to the environment. Hard boundary conditions mean that if a move has been selected which would shift an individual outside the lattice, the move is not realized and the individual stays put. Since it would be dicult to operate with so many parameters we decided on giving to some of them constant values. It has been, however, checked that setting them free does not change the results in a signi cant way. Thus, we have decided on taking LG = 5, q = 4, and x(0) = 0:3. This leaves only two free parameters – the number of loci, c, of the genetic string and the maximum age M for an individual. For the simulations we took a 100 × 100 square lattice. Averaging was done over 25 independent runs.
328
A. Pekalski / Physica A 252 (1998) 325–335
3. Results The quantities we shall record and discuss are the average tness hf(t)i =
N (t + 1) ; N (t)
(2)
average adaptation N (t)
a(t) =
1 X aj (t) ; N (t)
(3)
j=1
average age (relative to the maximum one) N (t)
w(t) =
1 X wj (t) ; MN (t)
(4)
j=1
average density of the population x(t) =
N (t) ; L×L
(5)
where L×L is the number of all sites on the lattice (carrying capacity of the system). At chosen time steps we shall record the distribution of the phenotypes, i.e., frequency of appearing of a phenotype with a given number of 1’s and the distribution of adaptation on a scale from 0 (not adapted at all) to 1 (perfectly adapted) with a 0.1 step. Finally, we de ne the tness function, t(t), which in a traditional way [8], is connected to fertility rate, i.e., it is the average number of osprings an individual had in its lifetime. Fig. 1 shows the time dependence of the average tness and adaptation. It is evident that the initially totally random system (adaptation ≈ 0.5, highly oscillating average tness) becomes more ordered in time. It adapts to the environment and the oscillations diminish. The amplitude of the oscillations increases with the number of loci in the genetic string (Fig. 2). This seems natural, since a population of individuals with longer phenotypes have more degrees of freedom. Somehow similar oscillations are observed for tness, t(t) (Fig. 3) a quantity not discussed in Ref. [4]. In this case the oscillations are more irregular but much less damped. Populations of individuals with short genetic strings (c = 10) also show uctuations of tness, which were not observed for the average tness de ned in Eq. (2). Individuals with shorter genetic strings adapt faster to the environment (see e.g. Ref. [3]), hence their tness is higher than for individuals with longer genetic string. Fig. 4 illustrates changes in the distribution of phenotypes in time. In agreement with the common assumption of biologists, it always has a Gaussian character. Since the initial population (“founders”) was random, the initial distribution is centered at the middle. In time, the population adapts to the environment, shifting the maximum of the distribution towards larger values of the phenotypes, i.e., genetic strings with more 1’s (good genes). The process has an asymptotic character, tending to minimize
A. Pekalski / Physica A 252 (1998) 325–335
329
Fig. 1. Time dependence of the average tness and adaptation. Length of the genetic string c = 20. Maximum age M = 5. Average over 25 independent runs.
Fig. 2. Time dependence of the average tness. The parameter is the number of loci. Maximum age M = 5. Average over 25 runs.
330
A. Pekalski / Physica A 252 (1998) 325–335
Fig. 3. Time dependence of the tness; c = 20, M = 5. Average over 25 runs.
Fig. 4. Distribution of phenotypes at chosen time steps; c = 20; M = 5.
A. Pekalski / Physica A 252 (1998) 325–335
331
Fig. 5. Time evolution of the maximum in the phenotypes distribution (from Fig. 4).
the number of defective loci. The distribution curves are in time getting more and more narrow. Note that this adaptation happens even though we did ignore here the possibility of mutations. It means that in the present model a population with initially random distribution of phenotypes can achieve perfect adaptation to the environment. This is reached only asymptotically, i.e., for an in nitely large population and after in nitely long time (see Fig. 5). In this model evolution is driven by two forces only – “survival of the ttest” rule and sexual reproduction. No mutations are needed. One may also construct the time dependence of the distribution of adaptation (Fig. 6). Like in Fig. 4 the distribution has always a Gaussian character and is shifted with time towards higher values, i.e., there are more and more better adapted individuals in the population. This adds arguments to the biologists’ claim that the distributions encountered in nature are generally Gaussian. All that, together with the adaptation curve, Fig. 1, indicates that the population is undergoing a homogenization process forced by the selection of the environment. Maximum age is one of the factors in uencing the dynamics of population in our model. Fig. 7 shows that both the amplitude and the period of average tness oscillations depend on it. The period is simply equal to the maximum age. This is dierent from the observations of Doebeli [4], where the period depended only on the number of loci. One should, however, bear in mind that the two models are not identical.
332
A. Pekalski / Physica A 252 (1998) 325–335
Fig. 6. Distribution of average adaptation at chosen time steps; c = 20, M = 5.
Fig. 7. Time dependence of the average tness parametrized by the maximum age M . Number of loci c = 20.
A. Pekalski / Physica A 252 (1998) 325–335
333
Fig. 8. Time dependence of the concentration, relative average age of an individual and average tness; c = 20, M = 5; Single run.
To better understand the dynamics of the present model, let us inspect one run on a restricted time scale, e.g. from 4 to 20 MCS for M = 5 and c = 20, Fig. 8. The evolution of the population is an interplay among two “killing” mechanisms – bad adaptation and reaching the maximum age – and one procreation mechanism of putting osprings on available empty sites within a limited area. Let us start the process with a large population (high concentration) at time step 8 in Fig. 8. The number of available sites for osprings is in this case very small, hence the number of produced osprings is small too. One may say that the survivial rate for a new born is very low. As a result, in the next time step the concentration of individuals falls down. Since few osprings are born, the average age is mounting and reaches maximum (time step 9). In the following step there will be more individuals eliminated because of their age, which leads to further diminishing of the concentration (time step 10). There are, however, more empty sites, liberated in the previous step, hence more osprings are born, which leads to lowering of the average age and to the increase of the concentration (time step 11). Individuals are now mostly eliminated by bad adaptation. There are still empty sites, hence the concentration grows (time step 12). The number of empty sites diminishes, hence the number of osprings is reduced, average age increases and the concentration reaches maximum (time step 13) and the process returns to the starting point, with however dierent values of the parameters.
334
A. Pekalski / Physica A 252 (1998) 325–335
Fig. 9. Comparison of time dependence of the average tness in two models – with mobile and immobile individuals (see text); c = 50, M = 5.
So far our model, although similar, did not yield the characteristic “saw-like” oscillations of the average tness, found by Doebeli [4]. Let us now modify the present model by dropping the requirements that in order to breed, an individual has to move and that the partner is chosen from nearest-neighboring sites only. In this way we reduce the role played by inbreeding and make the model more mean- eld-like. The results for the average tness dependence on time, together with the previous ones, are shown in Fig. 9. It is evident that the oscillations are much stronger and less damped. The oscillations have indeed some resemblance to the ones obtained by Doebeli [4]. The dierence is that here there is still a weak damping and the period of the oscillation does not depend on the number of loci, but only on the maximum age. The crucial point in passing from strongly to weakly damped oscillations is neglecting the need for a parent to move. It seems, therefore, that a diusive motion of individuals permits a better use of empty sites, hence it reduces the oscillations. 4. Conclusions We have presented a model of population dynamics based on Monte Carlo simulations. The model is similar, although not identical, to a mean- eld-like one considered by Doebeli [4]. The global characteristics like average tness or average adaptation,
A. Pekalski / Physica A 252 (1998) 325–335
335
etc., are determined in our model by the behavior of individuals, i.e., mobility and the selection rules. The tness function was de ned as an average number of osprings produced by an average member of the population. The population is under a selective force from the environment. We have found, like in Ref. [4], the oscillatory character of the time dependence of the average tness. The oscillations are damped, which indicates that the selection drives the population towards a homogeneous one. This seems quite natural since there are no mutations in our model. The damping and the period of oscillations depend on the maximum age allowed for an individual. The distribution of phenotypes is always Gaussian. The same is true for the distribution of adaptation of individuals to the environment. If we simplify the model by immobilizing the individuals and choosing partners at random from all population, not only from nearest neighbors, the oscillations of the average tness become much more pronounced and less damped. This is similar to the results produced in the mean- eld-like model in Ref. [4]. Even in our simpli ed model we did not nd the double-peak distribution of phenotypes. It was at each time step Gaussian. Acknowledgements I am very grateful to Dr. S. Cebrat for his valuable comments. References [1] [2] [3] [4] [5] [6] [7] [8]
P. Derrida, P. Higgs, J. Phys. A 24 (1991) L985. F. Manzo, L. Pelitti, J. Phys. A 27 (1994) 7079. I. Mroz, A. Pekalski, K. Sznajd-Weron, Phys. Rev. Lett. 76 (1996) 3025. M. Doebeli, Evolution 50 (1996) 532. J. Maynard-Smith, M. Slatkin, Ecology 54 (1973) 384. M. Slatkin, Ecology 64 (1980) 163. I. Mroz, A. Pekalski, submitted. T. Nagylaki, Introduction to Theoretical Population Genetics, Springer, Berlin, 1992.