International Journal of Stochastic AnalysisVolume 2013 (2013), Article ID 826321, 23 pageshttp://dx.doi.org/10.1155/2013/826321

Review Article

Simulating the Emergence of Mutations and Their Subsequent Evolution in an Age-Structured Stochastic Self-Regulating Process with Two Sexes

Charles J. Mode,¹ Candace K. Sleeman,² and Towfique Raj³

¹Department of Mathematics, Drexel University, Philadelphia, PA 19104, USA²Nokia Corporation, 5 Wayside Road, Burlington, MA 01803, USA³Division of Genetics, Harvard Medical School, Brigham and Women’s Hospital, Boston, MA 02115, USA

Received 28 August 2012; Revised 30 November 2012; Accepted 3 December 2012

Academic Editor: Peter Olofsson

Copyright © 2013 Charles J. Mode et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The stochastic process under consideration is intended to be not only part of the working paradigm of evolutionary and population genetics but also that of applied probability and stochastic processes with an emphasis on computer intensive methods. In particular, the process is an age-structured self-regulating multitype branching process with a genetic component consisting of an autosomal locus with two alleles for females and males. It is within this simple context that mutation will be quantified in terms of probabilities that a given allele mutates to the other per meiosis. But, unlike many models that are currently being used in mathematical population genetics, in which natural selection is often characterized in terms of parameters called fitness by genotype or phenotype, in this paper the parameterization of submodules of the model provides a framework for characterizing natural selection in terms of some of its components. One of these modules consists of reproductive success that is quantified in terms of the total expected number of offspring a female contributes to the population throughout her fertile years. Another component consists of survival probabilities that characterize an individual’s ability to compete for limited environmental resources. A third module consists of a parametric function that expresses the probabilities of survival in a birth cohort of individuals by age for both females and males. A forth module of the model as an acceptance matrix of conditional probabilities such female may show a preference for the genotype or phenotype as her male sexual partner. It is assumed that any force of natural selection acts at the level of the three genotypes under consideration for each sex. By assigning values of the parameters in each of the modules under consideration, it is possible to conduct Monte Carlo simulation experiments designed to study the effects of each component of selection separately or in any combination on a population evolving from a given initial population over some specified period of time.

1. Introduction

Comparatively little attention has been given to models on the evolution of age-structured populations in the extensive literature on evolutionary and population genetics. Much of the work on this class of models has, however, been reviewed and extended in the book by Charlesworth [1] on evolution in age structured populations. This book also contains a long list of references citing papers that are quite well known by, among others, Norton, Haldane, and Fisher, but it is beyond the scope of this paper to provide an overview of the content of these early works. It is of interest to note, however, that Norton began work on models of age structured populations as early as 1910, soon after Mendel’s laws of inheritance had become widely known. For the most part, models on the evolution of age structured populations that have been presented in this literature belong to the deterministic paradigm, but in this paper a special class stochastic models on the evolution of age structured populations will be the focus of attention.

This special class of stochastic processes is rooted in ideas from branching processes that also have an extensive literature dating back to about 100 years. Rather than attempting to cite papers from this extensive literature, books on branching processes will be cited which do contain extensive lists of literature. During the second half of the 20th century, a book by Harris [2] provided a stimulus for other workers enter the field of branching processes. In the period that followed the publication of this book, other books on this subject were published. Included in these books are those of Athreya and Ney [3], Jagers [4], and Mode [5]. In particular, the class of branching process under consideration in this paper is an extension of the general class of branching processes presented in Jagers, that contains a one type formulation of this process and Mode, which contains a multitype version of the general branching process. The multitype case of the general processes is the most useful when formulating processes applicable to evolutionary and population genetics, because this case accommodates not only mutations among the types but also differences in reproductive success among the types as well as the age structure of an evolving population. Other more recent books on branching processes are Asmussen and Hering [6], Kimmel and Axelrod [7], and Haccou et al. [8]. Each of these books has an extensive list of references that an interested reader may wish to explore.

The working paradigm underlying the class of branching processes to be formulated and studied in this paper differs markedly from that used in the books cited above. In these books, as well as the literature on branching processes published elsewhere, emphasis was placed on finding threshold conditions under which the population either became extinct or its total number, total population size, would increase without bound. Perhaps the most significant departure from the working paradigm of classical branching processes considered in this paper is that the process is self-regulating in the sense that its total size is dependent on the carrying capacity of the environment expressed in terms of parameters of survival functions, in which the parameter values may differ among genotypes. Because the survival of individuals during each time period under consideration is a non-linear function of total population size, the formulation being considered has nonlinear properties which arise in connection with the interaction of individuals in the population, as they compete for food and other resources. There are also other non-linear properties of the formulation under consideration that arise in connection with females selecting male sexual partners by genotypes that leads to the production of offspring that are necessary for the long-term survival of a population. Given a formulation that accommodates selection of mates by genotype, it is possible to study the effects of sexual selection on the long-term evolution of an age-structured population. For the most part, interactions among individuals were not accommodated in the formulations of classical branching processes in the books cited above, but the more recent book by Haccou et al. [8], where a one type a density dependent process is formulated in terms of, among other things, the well-known logistic function that played a role in the development of chaos theory that may be of interest to readers, but it is beyond the scope of this paper to present further details here. It should also be mentioned that, unlike the classical analysis of branching processes that often focuses on limit theorems, it is possible to study the long transient period of evolution a process before the process converges to some limit by using Monte Carlo simulation methods.

There is another aspect of the formulation to be considered in this paper which is also a significant departure from the working paradigm of classical branching processes, namely, the embedding of systems of nonlinear difference equations in the stochastic processes, which will be discussed in detail in a subsequent section and will be referred to as the embedded deterministic model. After software had been written to provide a computer implementation of the embedded deterministic model, given a set of numerical assignments to the parameters, it was possible to conduct computer simulation experiments, such that 10,000 to 20,000 years of evolution of an age structured population may be simulated with a computer running time from one to four minutes. Given this rapid computer response time, it was easy to determine whether a particular set of assigned numerical values of the parameters belong to a set of points in the parameter space, such that the population becomes extinct or grows to a total size where environmental resources limit further growth. Consequently, in this paper the task of finding criteria that partition a multidimensional parameter space into regions such that, according to the embedded deterministic model, the population either becomes extinct or its total size approaches some limit determined by the environment will not be under taken, but it is hoped that some analysts may be stimulated to undertake this task. After interesting regions of the parameter space have been found while using the embedded deterministic model, a Monte Carlo simulation experiment is conducted in which sample realizations of the sample functions of the stochastic process are computed using the same parameter values as those used in an experiment with embedded deterministic model to test the extent to which the predictions of deterministic are consistent with those of the process. As will be shown by examples, the predictions of the embedded deterministic process are not always consistent with those of the process, which is another facet of how the working paradigm of this paper differs from that of classical branching processes. It should also be mentioned that the formulation of an age structured process is essentially a complete revision of the age structured process considered in chapter 12 of Mode and Sleeman [9], where all reported computer experiments were based on applications of the embedded deterministic model.

2. Genotypes, Gametes, Mutation, and Types of Matings in an Age Structured Population

Only three genotypes with respect to two alleles, and , at some autosomal locus will be considered, and let denote the set of three genotypes. To simplify the notation, elements of will be denoted by the symbol , and to distinguish the sexes, female and male genotypes will be denoted by and , respectively. In what follows, the genotypes in the set will sometimes be denoted by 1, 2, and and will be made clear what sex is under consideration. The number of age classes under consideration is for both sexes, and the ages of females will be denoted by , where denotes infants or newborns. Similarly, ages of males will be denoted by . Let denote the minimum age at which fertile females may bear children, and let denote the maximum age at which they fertile. Let denote the ages of males, such that they may sire offspring. From now on such males will be referred to as fertile. It will be assumed for the sake of simplicity that .

Mutations among the alleles and , which will be assumed to occur during meiosis in both sexes, will be formulated in terms of conditional probabilities in the matrix where is the conditional probability that allele mutates to allele per meiosis. The conditional probability that allele does not mutate is . Similarly, is the conditional probability that allele mutates to per meiosis and . It will be assumed that these conditional probabilities hold for both sexes as well as throughout the fertile ages for both females and males.

The next step in the formulation of the model is to derive a gametic distribution for each of the three genotypes, which will be symbolized as follows. Let denote allele , and let denote allele , then the set of two alleles under consideration may be denoted by . Similarly, the set of three genotypes will be denoted by . When there are mutations among the two alleles, each of these genotypes may produce a gamete containing an allele . Let denote the conditional probability that genotype produces the gamete . By way if an illustration, a formula for the conditional probability will be derived. It will be assumed that the allele on the left in the symbol was contributed by the maternal parent of an individual and the one on the right was contributed by the paternal parent. The conditional probability that the allele on the left in the genotype is contributed to the gamete pool of the population is . Similarly, the conditional probability that the allele on the right is contributed to the gamete pool is . By adding these probabilities it follows that . By a similar argument, it can be shown that . Furthermore, by continuing this line of reasoning, it can be shown that ,  ,  , and .

A mating between a female and male of genotypes and , respectively, will be denoted by the mating type . Because there are three genotypes of each sex, it follows that there are 9 types of matings. Throughout this paper the set of mating types will be denoted by in the order indicated. For every mating type , let denote the conditional probability that a mating of type produces an offspring of genotype . By assumption females and males have the same gametic distributions, and it also assumed that female and male gamete combine independently. From these assumptions, it follows that for all couple types and genotypes . In computer implementations of the model under consideration, these conditional probabilities are computed in the order indicated in the set .

As discussed in Mode and Sleeman [9] in chapter 12, a module for couple formation, depending on the ages of the female and male, is an integral part of the formulation of an age dependent stochastic population process, then the number of couple types may become very large, which leads to the necessity of processing large arrays in computer implementations of the process. The processing of large arrays in a computer, in turn, becomes problematic, unless some scheme is adopted to reduce the number of couple types. It is of interest, therefore, to formulate a two sex age dependent population process in which sexual contacts may occur between females and males but without partnerships or couples which may last for long time periods. Accordingly, the purpose of this section is to formulate a process that includes sexual contacts between females and males but not partnerships of females and males of long duration but are transient and may vary from year to year among age groups of fertile females and males. This approach mating system modelled by this approach may be justified by assuming this type of polygamy was part of the evolution of the human species before the idea of long-term marital partnerships evolved.

3. Births in a Two Sex Age Dependent Population Process without Couple Formation

When an age structured population is under consideration, individuals in a population will also be classified by type. If for example, a female of genotype is of age , then her type will be denoted by . Similarly, a male type will be denoted by , where is his age and and belong to the set of genotypes. At this point in defining the components of the model, it should be mentioned that age and time will expressed in terms of a lattice: , where is some time unit. For example, if the evolution of a human population is under consideration, then usually will be a year. Hence, from now on . For every , let denote 3 by matrix valued function with 3 rows and columns, where denotes the number of females to type at time in a population. The 3 by matrix is defined similarly for males. To reduce the size of the arrays to be processed in a computer, it will be assumed that sexual contacts between females and males do not depend on the age of the male but only on the frequency of the male’s genotype. For , let the random function denote the number of males of type in the population at time . Observe that the random function is the number of individuals of genotype and age in the males population at time . Then it follows that, the random function is the total number of fertile males of genotype in the population at time . Therefore, the frequency of fertile males of genotype in the population at time is given by the random function for , where and if , where .

With a view towards minimizing the size of arrays that need to be processed in a computer, it will be assumed that for each fertile female of age expresses preferences for males sexual partners by genotype. To simplify the notation, denote the three genotypes under consideration by ,  , and , and let denote the matrix of acceptance probabilities for females. For example, is the conditional probability that a female of genotype 1 finds a male of genotype 1 acceptable as a sexual partner, and in general, given a female of genotype ,   is the conditional probability that she finds a male of genotype acceptable as a sexual partner. For the sake of simplicity, it will be assumed that this matrix does not depend on the age of a fertile female. If for all pairs , then, by definition, females select male genotypes at random for sexual contacts. But, if, for example, is greater than the other two probabilities for all , then all females prefer males of genotype 3 as sexual partners. This is an example of what is called sexual selection.

Given these definitions and an application of Bayes’s rule, it follows that is the conditional probability that a female of genotype has a sexual contact with a male of genotype during the any time interval , and let denote a vector of these conditional probabilities. At time , for , let the random function denote the number of fertile females of type in the population at time , and let the random function denote the number of females of type who have a sexual contacts with a males of genotype during the time interval . Then, let denote a random vector whose elements are the indicated random functions. It will be assumed that this random vector has a conditional multinomial distribution with index and probability vector . In symbols, for and . More precisely, the components of the vector are Observe for each the elements in the vector for are the random number of females of genotype and age that have sexual contacts with a males of genotype . Altogether there are 9 types of sexual contacts when there are 3 genotypes of each sex as was shown in the previous section. The random numbers of these 9 types of contacts during any time interval is given by the random vector Sexual selection, whereby females of all genotypes prefer males of one genotype over the others, is thought to be one of the driving forces of natural selection. The process just described accommodates this component of natural selection. It should also be mentioned that this simple mating process just described not only simplifies the computer implementation of the two-sex model under consideration, but it may also be a plausible model for the evolution of a population prior to the time long-term monogamous sexual partnerships evolved.

Reproductive success is also thought to be a component of natural selection. For each fertile female, who has sexual contacts with a fertile male, this component will be characterized by the parameters denoting the expected number of offspring each fertile female type contributes to the population during a time interval for . Human females are most fertile during their twenties, and the value of starting with will increase with up to about age 25 and then decreases to some value at age , the age at which a female is no longer fertile. Let denote the expected number of offspring produced by a female of genotype throughout her fertile ages. Then, for every , When assigning values to the parameters of a model in the process of desiging a computer experiment, suppose one assigns a value to each for . For example, for the case of a population with high fertility and no selection among the three genotypes, then one could assign the value for . Whereas, by way of an illustrative example, if it is assumed that selection is such that increases with , then the assigned values of the parameters could be ,  , and , so that, from the point of view of natural selection, genotype 3 would have a selective advantage over the other two genotypes.

Given (15), a question that naturally arises is how should one distribute the assigned value over the fertile ages, such that (15) holds? One approach to answer to this question is to find a probability density , such that Then, if , it follows that for all , so that (15) is satisfied. The density function may be called the fertility distribution of females during the fertile ages.

Among the many choices of the density function is a truncated Poisson distribution with parameter . It is well known that the mode of a Poisson density function is at the value when is a positive integer. If, for example, a female becomes fertile at age , then, under the assumption that her peak in fertility is at age 25, it follows that one should choose , so that the mode of the distribution is age at age . Next suppose that the maximum age a female is fertile is . Then the number of numbers between and starting at and ending in 40 is . To create a truncated Poisson distribution, compute the numerical values of the density function for . The next is to compute the normalizing constant Let denote the probability density function of Poisson distribution truncated at age 40. Then, for for . Other choices of the density of the distribution of the age of child bearing could be made, but, for the sake of simplicity, in all the computer experiments reported in this paper, this distribution will be that of the truncated Poisson described above.

In what follows, a well-known property of the family of Poisson distributions will be used repeatedly. Let for be a collection of independent Poisson random variables with expectations for all . Furthermore, suppose each of these random variables takes values in the set of nonnegative integers . Then, the random variable has a Poisson distribution with expectation (parameter) . It should be remarked that this property is easy to prove using probability generating functions. In what follows, it will be assumed that the number of females of age with sexual contacts of type produce offspring independently, and, moreover, females of different ages also produce offspring independently. It should also be mentioned that all conditional distributions and expectations are conditioned on the state of at time , so that the idea of conditional independence is in force.

In principle, under the assumption that the number of offspring produced by an individual female of genotype follows a Poisson distribution with expectation , then given a realization of the random variable and the expected value , it follows that the random number of offspring produced by females of genotype and age that had contacts with a males of genotype would follow a Poisson distribution with expectation . Then, because it is assumed that females of different ages also produce offspring independently, it follows that the random denoting the total number of offspring produced by females of genotype having male sexual partners of genotype follows a Poisson distribution with parameter during the time interval for all nine combinations of sexual partnerships .

To simplify the notation, denote the sum in parentheses by . Then has the simpler form From this formula, it can be seen that if , then the expected number of total number of offspring produced by partnerships of type during the time interval would exceed , indicating that if this relationship held forward in time, then the expected number of offspring produced by partnerships of type would increase over time. On the other hand, if and if this relationship held forward in time, then eventually mating types of the form would disappear from the population. It should also be noted that if there are no individuals of genotype in the population at time , then for all . It is of interest to observe that if the initial population consists only of individuals of genotype for both sexes, then initially for , until a sufficient numbers of mutant genotypes and have arisen as a result of mutations before for and . Moreover, if the initial number of females of genotype is small, then population may become extinct before the mutant genotypes and appear in the population. It is also of interest to observe that if for some all individuals, the female and male population, are of genotype , then except for the case .

Therefore, when writing software to simulate realizations of Poisson random variables, zero valued parameters need to be accommodated. Symbolically, let denote a procedure for simulating a realization of a Poisson random variable with parameter , given the state of the population at time . Let the random variable denote the number of matings of a female of genotype with a male of genotype during the time interval . Then, if , then , but if , then indicating that the random variable has a conditional Poisson distribution with parameter . When the parameter is large, then it is well known that the distribution of may be approximated by a normal distribution with expectation and variance .

In general, consider a random variable with a Poisson distribution with parameter , and suppose is large. Then, where is a standard normal random variable with expectation and variance and means approximately in distribution. When , the sum on the right is always positive, but if , then this sum may be negative, and, therefore, may be negative. But, every realization of a Poisson random variable must take a value in the set nonnegative integers. A question that naturally arises is if we set , where is the greatest integer function, then for what values of will the sum with sufficiently high probability? Observe that . Thus, , if and only if.

For any value of , it follows that because the standard normal distribution is symmetric about 0. Let be a trial value. Then, it can be shown by either computing or by consulting a reliable table for the standard normal distribution that . Thus, Therefore, for all such that , it follows that It seems reasonable, therefore, for all to use the central limit theorem approximation when simulating a realization of a Poisson random variable. When simulating realizations of many Poisson random variables, it is most efficient to use the normal approximation whenever it is reasonable, because each realization requires only one call to a procedure for simulating a realization of a standard normal random variable. Consequently, the procedure just described was implemented in the software, but, of course any investigator would be free to choose any trial value of that suits his goals.

For those readers who may be uncomfortable with this choice of the Poisson distribution for biological reasons, it would be straight forward if, for example, an investigator chooses to use negative binomial distribution for the offspring distribution in a branching process formulation. A procedure for simulating realizations of random variables would be more complicated using this two-parameter distribution than that described for the one parameter Poisson distribution. It could, however, be accomplished with a little more effort from some investigator by writing the software to accomplish the task in a programming language of his choosing. But, just as in the case of using a Poisson distribution, at some point in the development of the software, a central limit theorem would need to be invoked to simulate realizations of sums of conditionally independent random variables. If a reader is interested in a review of parametric offspring distributions that have been used extensively in demographic studies of human populations, it is suggested that the book Mode [10] be consulted.

Given a realization of the random variable , let the random vector denote the random number of offspring produced partnerships of type of each of the genotypes . Then, the vector has a conditional multinomial distribution with index, sample size, , and probability vector derived in the previous section. In symbols, for every mating type

Let the random variable denote the total number of offspring of genotype produced by females during the time interval for . Given that the probability that an offspring is a girl is , then and the number of females with genotype born during this time interval . Therefore, the number of boys of genotype born during this time interval is for .

In general, for and at time , let the random functions and denote, respectively, the number of females of type and the number of males of type in the population at time for . Let denote a vector with the components for for females, and define the vector similarly for males. These column vectors contain the number of females and males by genotype, respectively, born during the time interval . Then let denote a matrix with the column vectors for , and define the vector for males similarly. Then, the matrix for females at time is given by the matrix where the, indicates that the matrix is catenated onto the vector to construct a matrix. The matrix is constructed for males in the same way. Observe that at time females of age are in column of the matrix . Thus, in the construction of the matrix , females of age at time are dropped from the population. But, because very few females will survive to age , this loss will be negligible. This remark also applies to males. In the next section, a parametric model of a birth cohort survival function, governing the survival of individuals in a birth cohort as a function of their age will be described.

The next component of the age structured process under consideration is that of taking mortality into account in terms of condition survival probabilities by age. Before moving on to an overview of the survival component of the age structured process under consideration, it is appropriate to mention that the book Mode [10] also contains a description of stochastic models of human reproduction expressed in terms monthly estrous cycle of human females along with such factors as abortions, spontaneous and induced age of marriage that other factors involving the control of the number of births experienced by a female during her fertile period. Even though the inclusion of these details was useful in the quantifying of the short-term effects of family planning programs to include them in age structured evolutionary process under consideration with a long-term evolutionary perspective seems inappropriate for model that is in an early stage of development.

4. Survival of Individuals in a Birth Cohort as a Function of Age

Another component of natural selection that can be accommodated in the evolution of the age structured population under consideration is that of mortality of individuals. Let denote the probability that an individual in a birth cohort, who by definition was age at birth, survives to an age of at least . In symbols, if is a random variable taking values in the set of nonnegative real numbers, where is the set of real numbers, denoting the life span of an individual, then This function has the property that and is a monotone decreasing function, as increases and in the limit The purpose of this section is to derive a formulas for a survival function in terms of parametric risk functions.

For many species of animals, offspring are at high risk of death following hatching or birth, but as age increases the risk of death decreases. Therefore, the risk function for deaths of individuals for whom this risk function applies will be assumed to have a latent risk function of the exponential form as follows: where , and and are positive parameters. The integral of this latent risk function is for , and, by using well known methods for expressing a survival function in terns of an integral of the risk function, it follows that the survival function corresponding to this risk function is In what follows, it will be assumed that this survival function is in force for ages , because after age 30 the risk of death function will be assumed to be increasing. From now on all parameters of survival functions will depend on genotypes and sex of individuals, but simply the notation the genotype and sex of individuals will be suppressed except in cases when clarity is an issue.

Observe that , so that as it should. The integral will need to be modified to accommodate the conditions that . Let and let the modification of be defined as Let denote the modified survival function. Then it follows that , so that , which is the probability that an individual born at age survives to at least age 30. Thus, if one has some prior knowledge of the probability , then . Given this estimate of , a task that remains is that of finding a value for the parameter . Next observe that is the conditional probability that an infant born at survives to age year. Therefore, if one has some prior knowledge of the probability , by knowing the parameter could, in principle, be determined by finding a value of that would satisfy the equation . This is a nonlinear equation in the unknown parameter , and it may be possible by plotting as a function of to find a value , such that .

But, there is at least one other approach to find a value of . Observe that is the distribution function of a exponentially distributed random variable with expectation Thus, the parameter will determine the speed at which deaths occur. Small values of will correspond to longer infant survival times, while large values of will correspond to shorter survival times. Such ideas may be quantified by assigning trial values to to find an estimate of the form . For example, a plausible value of this expectation may be chosen as , which yields as a preliminary estimate of . The parameter may also depend on the sex and genotype of an individual and will be denoted by the symbols and for females and males of genotypes and , respectively. Given a value of , one would have to check whether is a plausible estimate of , the conditional probability that an infant born at age survives to age year. If the estimate of is not plausible, too low for example, then it would be necessary to try another estimate of .

To accommodate accidents that may occur throughout the life span of an individual, the latent risk function for this component of a survival function will have the simple form for all , where is a positive constant. In this case, the integral of the risk function is for . Therefore, for , the latent survival function has the simple form A useful trial value of , based on period studies of human mortality, is about 0.001. This component of the model is often referred to as the Makeham component and is named after a 19-century investigator that introduced this component.

A third two-parameter latent risk function, due to Gompertz (19th century), deals with risks of deaths at the older ages. Let and be positive parameters. Then, it will be assumed that for the latent risk function has the form It will be assumed that this risk function is in force for ages . Observe that, as it should, this risk function increases as age of an individual increases, and, by assumption, the risk of death increases exponentially with increasing age. The integral of this risk function has the form for . Therefore, the latent survival function for this component is for .

By applying a general but standard formula that a density is the risk function times the survival function, it can be seen that the probability density function of the Gompertz distribution has the form for . Although this distribution may be derived from intuitively appealing assumptions, it is more difficult to handle from a mathematical point of view than some other distributions that arise in probability and statistics. Nevertheless, because many advanced mathematical functions are now available to research workers in such software packages as MAPLE, MATHEMATICA., and MATLAB, an outline of the mathematics used in analyzing the Gompertz distribution seems appropriate.

Because the parameters and do not have obvious statistical interpretations, such as an expectation or variance, it is difficult to assign tentative values to them. Quite often, however, there is some feeling about the modal age of death for those who survive to old age. Let denote the mode of the Gompertz distribution. Then by using elementary calculus to find the maximum of the density , it can be shown that the equation formalizes a connection among the parameters ,  , and . In particular, if is assigned a value and is known, then is determined. But, to find a plausible value of , more input is needed.

Let denote a random variable with a Gompertz distribution. Then, after considerable analysis, it can be shown that the exact formula for the expectation is where is Euler’s constant. Furthermore, the exact formula for the second moment is

When is small, then and the above infinite series may be neglected. Thus, the approximations hold for small . Therefore, a formula for the approximate variance of the Gompertz distribution is Equivalently,

It will be instructive to present a simple numerical example, illustrating the use of these approximations. Suppose, for example, one wished to construct a Gompertz distribution with a mode years, and suppose the standard deviation of this distribution is . Then, this formula yields the approximation . Furthermore, by using the formula connecting ,  , and , it can be seen that this formula yields the approximation . Moreover, Because of the approximation of the variance is , it appears that the parameter will belong to the interval for plausible values of , so that it can be seen that when the mode of the distribution is sufficiently large, the parameter will be small. As can be seen form the above discussion, when is small, the approximation for will yield good results. A more extensive account of parametric forms of risk functions for death as well as historical references has be given in Mode and Sleeman [11, Section  13.2], but the details will be omitted here.

In the class of self regulating branch process under consideration, the effects on the mortality of individuals are formulated as a function of total population size at time . The formula for computing a realization of the random function will depend on the class of stochastic process under consideration. For the class of branching processes under consideration, as stated in a previous section, let denote the number of females of type in the population at time , and let the random function be defined similarly for males. Then, let the random function denote the total number of females in the population at time , and let denote the corresponding random function for males. Then, is total population size at time .

It will be assumed that the probability for all members of the population at time survive to time is governed by a Weibull type survival function depending on total population size of the form At this point in the discussion, it is appropriate to discuss the rational for choosing values of the parameter . For example, suppose the carrying capacity of the environment is individuals. Then, may be chosen as . As indicated with regard to all parameters under consideration, the positive parameters and will depend on sex and genotype. It follows, therefore, if one genotype has a smaller value of the parameter , it will be able to survive in larger populations and will thus have a selective advantage over the others. From now on the function will play a role similar to that of an integral of a risk function as illustrated above.

But, just as the total fertility of a female was distributed over her reproductive ages, a density function will needed to distribute over the ages of individuals in the population at time , which gives rise to the question as to how such a distribution may be chosen. A plausible answer to this question is that it may be constructed from the risks functions as functions of age as illustrated by the parametric risk functions introduced above. Observe that the larger the values of a total risk functions of an individual at these ages, the greater will be the proportion assigned to the total contribution of the Weibull component.

Thus, if the Makeham component is included in the risk function for the ages , then the total risk function for these ages is For ages the total risk function is Let be defined by the sum Then, the density for distributing over the ages may be chosen as for . Given these definitions, the proportion, , for age of the total will be chosen as for .

The next step in the formulation of the mortality component of the process is to express the final survival function in terms of integrals of risk functions. Thus, for ages the sum of these integrals is Recall that is the integral of the constant Makeham risk function, namely, . For ages , Therefore, the survival function for all ages is

The last step in formulating the mortality component of the process is to derive a formula for the conditional probability that an individual that is alive at age survives to age . Let denote this probability. Then, From this formula, it follows that in order that , the condition must be satisfied for all . Equivalently, the function must be strictly monotone increasing as a function of . It seems likely that because all parameters will be small, the component for population density will be small in relation to the others in . It seems reasonable, therefore, that if an investigator in a preliminary experiment ignores the component for population density, then it is likely that strict monotonicity of this function will be preserved when the population density component is included. To complete of the definition of survival function for , by definition , so that .

When checking the plausibility of a conditional survival probability in a preliminary experiment, the probabilities and require special attention. For example, so that the probability on the right should be consistent with a plausible probability that an infant dies in his first year of life. Similarly, because the risk function is increasing for all ages , conditional probability should satisfy the condition , so that the condition should also be satisfied. If this condition holds when the component for population density is ignored, then it seems likely, for reasons stated above, that this condition will also hold when the population density factor is taken into account.

The component of natural selection that needs to be considered in formulating the mortality component of the stochastic process under consideration is that of describing a procedure for simulating the number of individuals classified by age, sex, and genotype that are alive at time and survive to time . For example, let denote the conditional probability that a female of genotype and age at time survives to age at time . At time as before let denote the number of females of , and let denote the number of these females who survive to age at time . Then, by assumption, has a conditional binomial distribution with index, sample size , and probability . In symbols, Similarly, if is the number of males of type alive at time , then for .

5. An Embedded Deterministic Model in an Age-Structured Stochastic Process

As demonstrated in recent publications, Mode and Sleeman [9], Mode et al. (2011), [12, 13], by deriving formulas for the conditional expectation of any random variable at time , given the evolution of the process prior to , it is possible to derive a set of recursive nonlinear difference equations, such that, given the initial values of the random functions at time , estimates of the sample functions of the process can be derived for all . In previous publications dating back at least a decade, this derived set of equations has been called the deterministic model embedded in the stochastic process under consideration. The purpose of this section is to derive formulas for a deterministic model embedded in the age structured stochastic process under consideration. As will be seen in what follows, the derivation of formulas for the embedded deterministic model is closely related to procedures for computing Monte Carlo realizations of the process.

To fix ideas, suppose at time the elements of the random matrices and , denoting, respectively, the assigned numbers of females and males classified by age and genotype. All these initial numbers are nonnegative integers. Next consider the matrix valued conditional expectation Because is known, is known. Thus, is, as is well known that , the best estimate of in the sense of minimum mean square error. This procedure may be continued recursively by considering the conditional matrix valued conditional expectation But, this conditional matrix valued expectation is a random valuable, because is a random variable. However, if one has the estimate that has been calculated recursively, then is a reasonable estimate of the matrix valued random function for . This idea is the essential concept underling the procedure for embedding a recursive deterministic system in a stochastic process whose sample functions take values in matrices of non-negative integers. Deterministic estimates of the matrix valued random functions for males may be derived by an identical procedure for . Observe that the stochastic process and the embedded deterministic equations have the same parameter space, which is usually multidimensional.

Given these definitions, suppose that all the random matrix valued functions for males in (5), (6), and (7) are all replaced by estimates of the form starting in (5). Then, the estimate of the vector of the conditional probabilities in (9) have the form for . Observe that, because is constant for all nine pairs of matings types, the parameters do not need to be estimated. Thus, the estimate of the probability vector in (10) is

The next step in deriving the embedded deterministic model is to describe a procedure for estimating the random functions , which by definition is the number of matings of females of age and genotype with males of genotype , during the time interval . According to (11), the random vector has a conditional multinomial distribution index, sample size, , and probability vector in (10). As is well known, given the state of the population at time , the conditional expectation of the vector is the vector Therefore, let denote the estimate of the random function . Then, is an estimate of the random function for all pairs of mating types of females of age , where and are estimates, respectively, of the random functions and at time . From these definitions it follows that is an estimate of the random function at time , see (21). Let where is the density of the Poisson distribution defined in (20). Then, at time is an estimate of for all mating pairs of type , where is the expected potential number of total offspring produced by a female of genotype throughout her life span.

Let the random function denote the number of offspring produced by matings of type during the time interval . When it is not zero, then according to (23) this random function has a conditional Poisson distribution with parameter and sample size 1. Therefore, the estimate of this random function is The next step in deriving formulas for the embedded deterministic model is to estimate the random functions defined in the vector (28). Because, by assumption (29), the elements of this vector have a conditional multinomial distribution with index and probability vector derived in a previous section, it follows that an estimate of the random function is for all mating types and genotypes . Observe that the probability depends only on the probabilities of mutation, which, by assumption, are constant and, therefore, there is no need to estimate them during any time interval. From these results, it follows that is an estimate of the total number of offspring of genotype produced by females during the time interval .

Let denote the conditional probability that an offspring is a female. Then, according to (31) an estimate of , the number of females of genotype born during the time interval , is From this result it follows that is an estimate of the number of males if genotype born during the time interval . The estimate for the matrix at time is constructed by using the formula where is a vector with elements for . If a reader is interested in more details, the discussion pertaining to (33) should be consulted, because the procedure for constructing estimates of the random matrix involves replacing a call for random numbers by a conditional expectation. The construction of an estimate of the random matrix for males is identical of that for females.

The last step in constructing estimates of the state of the population at time , given estimates at time , is to replace the calls for random variables in (70) and (71) with conditional expectations of random variables with binomial distributions for those individuals who survive from to . Thus, for all genotypes and ages .

6. Overview of Procedures and Issues Involved in Designing Computer Experiments

Among the issues that arise as one begins the process of designing a computer simulation experiment is that of choosing a numerical value for each of the parameters in the multidimensional parameter space of the age structured stochastic population process under consideration. As can be seen from Section 2, the minimal age at which a female becomes fertile and may bear children as well as the maximum age that she is fertile must be assigned numbers. For most human populations, the age at which females become fertile varies around 15 years. In all computer experiments reported in this paper, however, variation among the ages when a female becomes fertile will be ignored, and for the sake of simplicity that parameter will be assigned the value . For similar reasons, the maximum fertile age for females will be assigned the value years. Another parameter that must be assigned a numerical value is , the greatest age an individual, female or male, may attain. For most of the computer experiments reported in this paper, this parameter was assigned the value years. Part of the rationale for choosing this value was the thought that in ancient populations the life span of individuals was short when compared with those of some modern populations, so that age 50 years would be a plausible maximum life span for an ancient human population. The other issue underlying this choice was the time required to complete a computer experiment covering an evolutionary time span of thousands of years. Simply put, the larger the value chosen for the parameter , the greater the time span required to complete a Monte Carlo simulation experiment consisting of some number of replications. It was found in preliminary experiments with the assigned value could be completed in reasonable lengths of time, given personal computer resources that were available to conduct the experiments.

The evolutionary driving force of mutation will play a basic role in all computer experiments reported in this paper, and, as can be seen from (2), the probability that allele mutates to allele per meiosis must be assigned a numerical value. Similarly, the probability that allele mutates to allele must be assigned a numerical value. The assigned values for these mutations will differ among the experiments reported in this paper; thus, in what follows, the values used in the experiments reported in subsequent sections will be sated in the initial computer input for each experiment.

In Section 3 on simulating the number of births that occur in a population in any year, it was shown that the birth process depends on two sets of parameters. The first set of these parameters is set forth equation (8) where a matrix of acceptance probabilities displayed that characterizes preferences of females of genotype for males of genotype as sexual partners for all pairs of sexual partnerships. The second set of parameters is , where the expected total number of offspring produced by a female of genotype during her fertile period from ages up to and including age . The values chosen for these two sets of parameters will also vary among experiments, so that no specific values will be given here. In all experiments reported in this paper, it was assumed that the probability that a baby was a girl was assigned the value , and by assumption, the probability a baby was a boy was assigned the value .

In Section 4, which is devoted to the survival of individuals in a population by age and from year to year, the two-parameter Weibull distribution was chosen to provide a framework to model the regulationof total population size. As may be seen from (61), a function of the form , where is total population size in year , and and are positive parameters that may differ by sex and genotype of an individual, must be assigned numerical values. Because two sexes and three genotypes are under consideration, the total number of alpha and beta parameters is 12. For the sake of simplicity, in all experiments reported in this paper, the values of all alpha parameters were chosen as , which was equivalent to assuming that all survival functions had a simple exponential form. This simplification of the formulation reduced the dimension of the parameter space to 6, but in order to provide a framework for simulating the evolution of a population that was, in part, dependent on the carrying capacity of the environment, it would still be necessary to assign to six beta parameters to accommodate the notion that one genotype may have a selective advantage over others by being a superior competitor for resources and was able to survive and reproduce in higher population densities. To take into account that individuals of genotypes may be more competitive users of environmental resources, and in all computer experiments reported in this paper, the values of the beta parameters will be given. Even though the survival capabilities may differ by sex, it was assumed, for the sake of simplicity, in all experiments reported in this paper that the assigned value of each parameter in the parameter space was equal for males and females.

In the component of survival function described in the first paragraphs of Section 4 for ages , the parameter , which may differ among the three genotypes and the two sexes, was among the parameters that must be assigned numerical values when designing a computer simulation experiment. As indicated in Section 4, this parameter has the interpretation that in a birth cohort for any year may be interpreted as potential probability that an individual in this cohort would survive to age 30. In all computer experiments reported in the experiments reported in this paper, this parameter was assigned the value for each sex and genotype. For the sake of simplicity, it was also assumed that the parameter had the same value, , for both sexes and genotypes within sexes. Thus, by assumption, in the computer experiments reported in this paper, no genotype or sex had a selective advantage over others with respect to the component of natural selection expressed in terms of this survival component of the formulation for ages .

Another aspect of the formulation governing the survival of individuals was that of the one-parameter Makeham component, taking into account deaths due to accidents at all ages. The positive parameter of this component was denoted by , and for each sex and genotypes within sexes, this parameter was assigned numerical value , which has also been used by others in connection with parametric models of human mortality. The technical details concerning this component are expressed in (42) and (43). The survival component in the formulation for ages was described in terms of the two parameter risk function in (44). As was shown in the discussion following this equation, the positive parameters and of the Gompertz distribution may be expressed in terms of the mode and standard deviation of this distribution. For each genotype within sexes, the mode and standard deviation of this distribution were assigned the values and years, respectively. In summary, the procedure just described, entailing the assignment of values of the parameters that were the same for both sexes and genotypes within each sex, was the implementation of an approach to quantifying the idea that forces governing the evolution of the age structured population were neutral with respect to survival by age.

The last component of the formulation that requires assignments of values, expressed in terms of non-negative integers, is the elements of the 3 by matrices and , denoting the number of individuals by genotype and age in the initial population at time . In all experiments reported in this paper, it was assumed that the initial population for both sexes consisted only of individuals of genotype , so that the genotypes and could arise in an evolving population only as a result of the allele mutating to allele , see (2). According to the “Out of Africa Hypothesis,” the people who now make up the current populations of Europe, Asia, Australia, and the Americas are descendents of a group of our species who emigrated out of Africa about 60,000 years ago. If a reader is interested in a detailed account of this hypothesis by a paleoanthropologist, it is suggested that the video course and book are Hawks [14] of consulted. It is also thought by many that this group of emigrants consisted of a rather small number of individuals with a total population size ranging from the hundreds to at most a few thousands. To simulate this initial founder population, whose size was uncertain, a random number generator was used to estimate the number of individuals of each age and sex for genotype 1, and the number 0 was assigned to all ages for individuals of genotypes 2 and 3. A broad objective of the computer simulation experiments reported in this paper is to study the evolution of the mutant genotypes and in an age structured population evolving from the ancestral genotype .

Various procedures will be used to statistically summarize data generated in Monte Carlo simulation experiments. These methods will not be described here, but if a reader is interested in a detailed overview of these methods, it is suggested that the paper of Mode and Gallop [15] is consulted. With a view towards complete openness as to the generation of random numbers used in the Monte Carlo simulation experiments reported in this paper, it is appropriate to mention the random number generator used in these experiments is that for 32-bit word computers designed by Deng and described in Mode and Gallop, see Section 2 of that paper. Further discussion of issues concerning random number generator may be found in the papers by Deng and Lin [16], L’Ecuyer and Touzin [17], and Deng [18]. But, unlike the 32-bit word computers that the Deng random number generator was designed for, the word length for computer used in the Monte Carlo simulation experiments reported in this paper was a 64-bits. Because the developers of the version 11 of the APL 2000 software used to write programs to implement the age structured process under consideration did not include a default 64-bit word random number generator, it was necessary for users to implement such a generator. However, to implement one of the 64-bit random number generators described in Deng would have required a rather long period of time to develop the necessary APL code in way, such that the times taken to complete Monte Carlo simulation experiments would be satisfactory. Consequently, this developmental task will be postponed to some future date. It is felt, however, that even if the software for a 64 bit random number generator would have been developed, from a qualitative point of view, the results of the Monte Carlo simulation experiments reported in this paper would not have been significantly different, due to, among other things, the long period of the random number generator used in these experiments, and it verified randomness properties.

Collections of human mortality data are also sometimes useful in choosing numerical input to age structured stochastic processes. In this connection, the collection of data in the book by Alderson [19] may be useful. A search for historical human mortality data on the internet may also yield informative information.

7. On the Emergence of Rare Mutations-Differing Deterministic and Stochastic Predictions

The objective of the computer simulation experiment reported in this section was to compare the predictions of the embedded deterministic model with those of the stochastic process with respect to the emergence of rare mutations. To simulate the occurrence of rare mutations, it was assumed that the conditional probability that the ancestral allele mutated to allele per meiosis was , and the conditional probability of the back mutation was assigned the value . These assignments were based on anecdotal evidence gleaned in conversations with microbiologists, who studied the evolution of a bacterial population for thousands of generations. The mutation studied in the experiment were rare in the sense that in previous experiments it was assumed that probabilities of mutation were in the interval per cell division or, in the case of diploids, per meiosis. It was also assumed that genotypes of individuals carrying the mutant allele had a reproductive advantage over individuals with the ancestral genotype. To quantify this idea, the potential number of offspring of females of each of the three genotypes throughout their fertile years were assigned the values . Given these assignments, it was expected that the population would grow rapidly so as to increase the likelihood of a rare mutation being realized in a Monte Carlo simulation experiment.

It was also assumed that mating was random in the sense that females of any genotype when choosing a male sexual partner had no preferences with regard to the genotype of her male partner. To quantify this assumption it was assumed that all the elements on the matrix of acceptance probabilities in (8) were assigned value 1, but it may be noted that any constant probability would suffice for the case that female mating preferences were neutral with respect to the genotype of the male. With regard to the carrying capacity of the environment, it was assumed that the beta parameters for each genotype within each sex had the constant value . It is of interest to note that this value is greater than the probability of a rare mutation, which raised the question as to whether the carrying capacity of the environment was sufficiently large to assure that at least one rare mutation would actually be observed in a Monte Carlo simulation experiment.

Given these assigned values of the parameters of the model, a preliminary experiment was conducted using the embedded deterministic model to simulate the evolution of an age structured population for 5,000 years. The numerical operations entailed in an experiment with the embedded deterministic model consist essentially of multiplications and additions, so that even for mutations with small probabilities of occurrence, the number of mutant genotypes will eventually exceed one for the case of high fertility considered in the experiment. As presented in Figure 1 there are the trajectories of total population size for each of the three genotypes in the female population for 5,000 years of evolution.

Figure 1: Trajectories of the total population size of the three female genotypes as computed using the embedded deterministic model.

As can be seen form this figure, the total size of the initial ancestral population of females, consisting only of individuals of ancestral genotype , reached the carrying capacity of the environment within 500 years of evolution. But, after about 3,000 years of evolution, the total size of population of individuals of the ancestral genotype began a steep decline as the total numbers of the females of the mutant genotypes and grew to significant numbers. Sometime before 4,500 years of evolution, the total number of individuals of genotype 2 starts a decline, as the total number of individuals of genotypes 3 rose to predominance in the population. This rise was due to the assumption that females of this genotype had a selective advantage over the other genotypes due to a higher level of reproductive success as quantified in terms of the numerical values of the parameters . That the trajectory for the total number of individuals of genotype 3 that was still rising near the carrying capacity of the environment suggests that this limit would be reached shortly after 5,000 years of evolution.

To check whether the predictions of a sample of the process would be such that the trajectories computed using the embedded deterministic model could be interpreted as measures of central tendency for the sample functions of the stochastic processes, a Monte Carlo simulation experiment was run with an evolutionary time span of 2,000 years with 50 replications. In preliminary experiments with 5 replications, it had been observed that within 2,000 years of evolution genotypes carrying rare mutant alleles had been realized in the simulated data. A decision was, therefore, made to run a confirmatory Monte Carlo simulation experiment with 50 replications of 2,000 years of evolution, which 28 took about 30 hours of computer time to complete the experiment. The reason for computing only 50 realizations of the process was purely practical, for it had been decided to compute 100 realizations of the process, and it would have taken about 60 hours of computing time to complete the experiment. Moreover, the risk of power outage, resulting in the loss of simulated data, during such a time period was significant for the home computer used to carry out the experiments reported in this paper. If a reader is interested in a more in depth discussion of the issues that arise in choosing the number of replications for a Monte Carlo simulation experiment, the paper by Mode et al. [13] may be consulted.

Whenever a Monte Carlo simulation experiment is carried out based on an age structured framework with several genotypes, a plethora of simulated data is produced, which gives rise to problems of presenting the results of such experiments in an informative manner. In this section, attention will be focused on the evolution of the population of individuals of genotype in which each individual carries two copies of the rare mutant allele . As presented in Figure 2 there are the , and trajectories of the total number of individuals of genotype 3 for 2,000 years of evolution from a founder population consisting only of individuals of the ancestral genotype .

Figure 2: Estimated quantile trajectories of the total number of females of mutant genotype 3.

As can be seen from this figure, significant numbers of individuals of mutant genotype appeared between 1,000 and 1,200 years of evolution, but, as can be seen from Figure 1, significant numbers of this genotype were not realized in the population until somewhere between 3,000 and 4,000 years of evolution according to the predictions of the embedded deterministic model. Consequently, this is a clear cut case in which the predictions of the timing of the emergence of a mutant genotype differs significantly as can be seen by comparing the trajectories of the stochastic process and the embedded deterministic model. If a reader is interested in other examples in which the predictions of the stochastic process differ significantly from those of the embedded deterministic model, the papers of Mode et al. (2011) [12, 13] may be consulted.

When an age structured model is under consideration, most of the simulated data arising in a Monte Carlo simulation experiment is in the form of arrays in which the number of individuals of each age is recorded. Furthermore, when many replications of realizations of the stochastic process are computed, it becomes necessary to statistically summarize the data in an informative manner. In the experiment under consideration the mean number of individuals for each age group and its standard deviation were computed. For example, let denotetheestimated mean number of individuals of age in the female population in year of a Monte Carlo simulation experiment, and let denote the estimated standard deviation. Then for year , the graph of the estimates is plotted in the upper panel of Figure 3 for ages and females of genotype 3. Similarly, in the lower panel of this figure the estimates of the standard deviations by age for this genotype are plotted.

(a)

(b)

Figure 3: Graphs of the mean and standard deviation by age for females of genotype 3 in year 300.

From these two graphs, it can be seen that the graph of number of individuals by age declines rapidly as age increases, which is a signature of rapidly growing population with quite high levels of mortality. It is also of interest to note that the number of individuals of the mutant genotype was present in the population, with numbers ranging from the low single digits to 25, within at least 300 years of simulated evolution, which illustrates that if one replied solely on the observing the graph of the total number of females with this genotype, this earlier appearance of the mutant genotype would have been missed. From the lower panel of this figure, it can be seem from the graph of the standard deviation by age is quite irregular for those ages representing the young, but these fluctuations level off age increases. The magnitude of the values of the standard deviation by age are indicators of the variations among the 50 replications of the process, particularly for ages less than 10 years. As presented in Figure 4 there are the graphs for the mean number of individuals by age along with the corresponding standard deviation for year in the Monte Carlo simulation experiment for females of mutant genotype 3.

(a)

(b)

Figure 4: Graphs of mean and standard deviation by age for females of genotype 3 in year 2,000.

From these graphs, it can be seen in the upper panel of the figure that, in year 2,000 of the Monte Carlo simulation experiment, the population of infants, , and ages near that of infants had reached a population size of nearly , that is, a size approaching the carrying capacity of the environment. In the lower panel, it can be seen that the standard deviations for these early years of life were near the large value , which is indicative of a rather high level of variation among the 50 replications of Monte Carlo realizations of the process for these ages. As age increases, the numerical values depicted in the graphs of the mean and standard deviations for the number of individuals decrease rather rapidly, which again was a signature of high levels of mortality in the population across the ages. It should be noted, however, that there is a steep drop in the graphs for both the means and standard deviations at age 30, where it was assumed that mortality beyond this age mortality was governed by a Gompertz type risk function. There are methodological reasons for this drop. A reader may recall that the potential probability that an infant survives to age 30 was assigned number 0.99. It seems plausible, therefore, that this potential probability should be lowered. On the other hand, this drop also suggests that different values for the two parameters of the Gompertz distribution may also help level off the steep decline in these graphs at age 30. While developing the formulation of the age structured stochastic process described in this paper, the most difficult problem encountered in this initial period of model development was the choices made in constructing and choosing the parametric components of the model governing the survival or mortality of individuals. It should also be noted that for ancient populations of man considered in the experiments reported in this paper, there is little or no data on mortality of individuals apart from the fossil record that suggest through measurements of radioactive that maximum ages of those individuals whose skeletons survived was in the range from 40 to 50 years. In passing, it should be mentioned that it is doubtful whether the changes in mortality parameters just suggested would have changed the differences in estimated timings for the emergence of rare mutations as observed in the realized predictions of the embedded deterministic and the stochastic process used in the simulation experiments reported in this section.

8. Sexual Selection in Evolving Age-Structured Populations

Among the many questions addressed by Darwin in his seminal work on evolution was that of the impact of sexual selection on the evolution of two-sex-species such as man, other mammals, and birds. A historical sketch of these ideas is beyond the scope of this paper, but if a reader is interested in the subject, a search of the internet may yield useful results, or one could also consult text books on population and evolutionary genetics, which are available in university libraries as well as for sale on the internet. For example, the book on population genetics by Hartl and Clarke [20] devotes some space to sexual selection. Essentially all the literature on population genetics in which models of sexual selection are presented are confined to deterministic paradigm, but in the book by Mode and Sleeman [9] a formulation of a two sex self-regulating branching process is the subject of chapter 11 along with reported results of some Monte Carlo simulation experiments. A limitation of the model considered in that chapter was that the age structure of an evolving population was not included in the formulation. Also this chapter contains the results of a Monte Carlo simulation experiment designed to study the rise of a recessive mutant allele in an evolving population when the character or trait it expressed was subject to sexual selection. In that chapter, sexual selection was characterized in terms of matrices of acceptance probabilities regrading both females and males as in (8). In Monte Carlo simulation experiments, it was shown that when it was assumed that the mutant genotype was preferred as a sexual partner over the other two genotypes, it would eventually arise to predominance in a population, but its steep rise to predominance would not occur until the total population size of individuals of genotype had reached a threshold. In this experiment the embedded deterministic model was a fair predictor of the central tendencies of Monte Carlo realization of the sample functions of the process, but in a related experiment, Mode et al. [12], the embedded deterministic model predicted a rare mutation that had sexual selective advantage would eventually become predominant in an evolving population, but in the corresponding Monte Carlo simulation experiment, the mutant recessive genotype did not appear, so that the predictions of the deterministic and stochastic models were not consistent.

In the Monte Carlo experiment reported in this section, the number of years of evolution under consideration was 10,000. Among the reasons underlying this choice was that about 10,000 years ago, in geographical regions of the world, such as Egypt and Mesopotamia, man began to change from a hunter-gather means of living to a more sedentary form of life, the development of agriculture that resulted in a more dependable source of food. At the outset it was known that if 50 Monte Carlo replication of an experiment were computed, the span if time taken to complete the experiment would be too long, given the private computer facilities available to carry out the experiment. So a decision was made to compute only 10 replications of 10,000 years of simulated evolution, which took about 30 to 40 hours of computer time. Because it was observed in some experiments not reported in this section that a significant number of individuals may be alive at age 50, a decision was made to set the maximum age at 80 years, which required the computer to process lager arrays that in turn would increase the expected time needed to complete an experiment. To take into account that sexual selection was operative in the experiment, the matrix of acceptance probabilities for females was chosen as

According to this selection of numerical values, the lower case allele , which arose from the mutation , was dominant with respect to sexual selection. Interestingly, in preliminary experiment with the embedded deterministic model, it was observed that if the values of the matrix were chosen, such that the values in the first and second columns were assigned the constant 0.1, and the values in the third column were assigned the constant value 0.9, then individuals of the mutant genotype would eventually become predominate in the population, but if these values were used in a Monte Carlo simulation experiment, then the computer time required to complete an experiment would be too long. Thus, for this reason the values in the matrix were chosen, because the study of the emergence of a dominant mutation would also be of interest in its own right.

By assumption, selection was neutral with respect to reproductive success so that the potential expected number of offspring contributed to the population by each female through out her lifespan was assigned the number 4 for all three genotypes. Observe that under this assumption the pace of evolution was expected to be slower than in the experiment reported in Section 6. To assure that an individual with a mutant gene would appear in a simulated population within an acceptable number of years, the probability of the mutation per meiosis was assigned the value , and the probability of the back mutation per meiosis was assigned the value . The carrying capacity of the environment for each genotype for both sexes was assigned the value, expressed in terms of the beta parameters, and had the constant value . All parameters not mentioned in the above discussion were assigned the same values as those in the experiment reported in Section 6. Moreover, it was assumed that the initial population of females and males was all of genotype and assigned the same initial numbers as those used in the experiment reported in Section 7.

As displayed in Figure 5 there are the trajectories of the three genotypes, expressed in terns of the numbers of individuals per genotype as functions of time, for the first 1,000 years as computed using the embedded deterministic model.

Figure 5: Deterministic trajectories of total numbers for each of the three male genotype under sexual selection.

As can be seen from these trajectories, from an initial population of 2,693 males of genotype , within 500 years of simulated evolution the mutant genotype starts to become predominant in the population, and by 1,000 years the number of individuals of this genotype has risen to over but has not yet reached the maximum carrying capacity of the environment. As one would expect, the number of individuals of genotype occurs in significant numbers in the time interval from 500 to 1,000 years but is far below the number of individuals of genotype 3 throughout this time period. An inspection of the simulated data in the years preceding 10,000 years indicates that a balanced polymorphism had been reached in which the genotype was predominant with genotypes and present in smaller numbers. Unlike the cases in which the recessive genotype is favored by sexual selection presented in Mode and Sleeman [9] and Mode et al. [12], there is gradual rise but no steep rise to predominance in the population of individuals of genotype .

From now on, due to limitations on space, attention will be focused on the predominant genotype 3 with respect to total population size for males of this genotype. As presented in Figure 6 there are graphs of the    and   trajectories as estimated from the Monte Carlo realizations of the process as well as the deterministic trajectory () for the total number of individuals of mutant genotype for the first 200 years of 10,000 years of simulated evolution.

Figure 6: Estimated    and   trajectories of the process and the deterministic trajectory () for the total number males of genotype 3.

As can be seen from this graph, after about 140 years of evolution, the number of individuals of genotype had risen to significant numbers, but after 150 years of evolution, the variation among the 10 realizations of the process increased as can be seen from the distances among the     and trajectories at year 200. Interestingly, however, at year 200, there is a large distance between the deterministic trajectory () and the of the stochastic process. This observation demonstrates that even with the high probabilities of mutation used in the experiment under consideration, the embedded deterministic model tends to underestimate the number of mutations that actually arise among the simulated realizations of the process in the first 200 years of simulated evolution. But, as shown in the next figure, after a sufficient amount of evolutionary time, the four trajectories in Figure 6 tend to coalesce.

As can be seen from this figure on the scale of the vertical axis, the distance between the and trajectories is close up to about 2,500 years of evolution. But, in the years following 2,500, the four trajectories coalesce and remain near one another at 4,000 years at evolution at least on the scale of the vertical axis. In the years following 4,000, it was seen, by inspecting the statistical summarizations of the Monte Carlo realizations of the process as well as that for embedded deterministic model, that the pattern observed at 4,000 years of evolution continued to 10,000 years. For the class of stochastic processes under consideration, the trajectory of the embedded deterministic model was quite typical in the sense that in the early years of an experiment the predictions of the deterministic model may differ significantly from those of the process, but in the long term, with the exception of a few cases, the trajectories of the stochastic process and that of the deterministic model tend to coalesce. But, the waiting time to coalescence may indeed be very long when expressed in terms of human life spans.

Because we are dealing with an age structured process, it is appropriate to view graphs of a simulated age distribution at some selected year of the experiment. As presented in the upper panel of Figure 8 ther is a graph for males of the mutant genotype , at year 200 of the experiment, and the lower panel contains the graph of the standard deviations for each age as estimated from the simulated data (Figure 7).

Figure 7: Stochastic and deterministic trajectories for the total numbers of males of genotype 3 for years from 1,000 to 4,000 of evolution.

(a)

(b)

Figure 8: Estimated mean age distribution and standard deviations by age for the number of males of genotype 3 in year 200.

From the upper panel of the figure, it can be seen that by year 200 of the experiment the number of males of the mutant genotype 3 had risen to sufficiently high numbers to form a nearly monotone decreasing mean age distribution that one would expect in a population evolving with high levels of mortality and fertility as realized in a large number of births and deaths each year. But as can be seen from the lower panel for the graph of standard deviations by age, there was a considerable amount of variation by age around the mean age distribution.

It is also of interest to view a three-dimensional version of the evolving mean age distribution for males of mutant genotype . As presented in Figure 9 there is a three dimensional view of the evolving age distribution for males of genotype 3 for the first 200 years of evolution.

Figure 9: Three dimensional view of the evolution of the mean age distribution for males of genotype 3.

As can be seen from this figure, the -axis denotes age, the -axis denoted time in years, and the -axis denotes the number of individuals. Unlike the conventional displays of an evolving age distribution, in which the numbers are normalized such that they all are in the interval , in this figure the number of individuals for each age is depicted. From an inspection of this figure, it can be seen that for the first of about 140 years of evolution the age distribution is a flat plane, indicating that the number of individuals of mutant genotype 3 was either zero at each group or a number too small to depict in a graph. However, sometime after 140 years of evolution the number of individuals of mutant genotype 3 increases in the population, and by year 200 the age distribution is typical of a fast growing population with high levels of mortality as is also shown in the upper panel of Figure 4.

A question that naturally arises is what is the three-dimension form of the three dimension standard deviation surface corresponding to the graph presented in Figure 9? Presented in Figure 10 is a graph of this evolving surface.

Figure 10: Three dimensional view of the evolving standard deviation surface for males of genotype 3.

The - and -axes are the same as those in Figure 9, but on the -axis the standard deviations are expressed in terms of units of individuals. Observe that shape of the standard deviation surface is very similar to that for the age distribution in Figure 9, but the range of values for the in Figure 10 is from 0 to 200; whereas in Figure 9 the range of the axis is from 0 to 600. Thus, at its maximum, the axis in Figure 10 is about only one-third as that for that axis in Figure 9.

A question that naturally arises is to what extent would the above figures change, if 50 to 100 realizations of the process had computed rather than the small sample of ten? It seems likely, based on other Monte Carlo simulation experiments, that from a qualitative perspective the estimated trajectories, as estimated from a Monte Carlo simulation experiment, would not change significantly, but if a larger number of realizations of the process were computed, it is likely more outliers would appear in the simulated data. Thus, for example, the at 200 years may have risen to a larger number than that displayed in Figure 6. Moreover, the standard deviations by age in the lower panel of Figure 8 would very likely be lager, indicating a greater range of variation among the simulated realizations of the process.

Another question that arises is what is the mathematical basis underlying the apparent stochastic-balanced polymorphism that was observed in the latter years of 10,000 years of simulated evolution? It can be shown that the two-sex age structured process under consideration is a Markov chain in discrete time with a state space consisting of pairs if matrices of non-negative integers, where the rows represent genotypes and columns ages of individuals and and females and males, respectively. In the absence of immigration, the state is an absorbing state, indicating extinction of the population. Given that extinction that did not occur, it appears that the process had converged to a quasi-stationary distribution of the Markov chain. No attempt will be made here to prove this assertion, but if the reader is interested in a derivation of a quasi-stationary distribution for a Wright-Fisher process with two absorbing sates, chapters 5 and 6 of Mode and Sleeman [9] may be consulted.

9. Discussion, Further Reading, and Future Experiments

In the two experiments reported in this paper dealing with selection, only one component of selection was dealt with in each experiment. For example, in one experiment dealing with the component of reproductive success, if females of one genotypes contributed to average a greater expected number of offspring than females of other genotypes, this selective advantage was sufficient for individuals of this genotype to become predominant in the population in the long run. On the other hand, if individuals of one genotype or a dominant phenotype in males was preferred by females as sexual partners for the case of sexual selection, this preference was sufficient for that genotype or phenotype to become predominant in a population in the long run. In both these experiments, all other components of selection were neutral in the sense that for each component not under consideration each genotype was assigned the same parameter values. In an experiment not reported in this paper, it was also shown that if individuals of genotype for either sex had a competitive advantage over the other two with respect to competition for resources, then individuals of genotype 3 would eventually become predominant in an evolving population.

A component of selection that has so far not been studied using the formulation under consideration is that of different survival patterns by age as formulated in terms of the parametric birth cohort survival function described in a previous section. For example, if females of one of the genotypes were very good mothers in the sense that most of their children survived to reproduce, then individuals of that genotype would have a selective advantage over individuals with other genotypes. It would be of interest to study this component of selection, but such computer experiments will be deferred to future research. Experiments dealing with combinations of components of selection would also be of interest. There is at least one other class of experiments that would be of considerable interest, namely, a case in which a human population evolves during a phase of evolution in which life is maintained by hunting animals, fishing, and gathering plant foods. Under such conditions the carrying capacity of the environment would be less than that considered in the experiments reported in this paper, in which it was tacitly assumed that agriculture had advanced to the point at which the environment could support population sizes of the order . In a hunter-gather society, however, a plausible estimate of the carrying capacity of the environment may be one or two people per square mile. As an illustrative example, suppose that the environment could support two people, a female and male, per square mile, then a population of 1,000 females and males would require an area of 500 hundred square miles to support such a population. Under such conditions, one would expect that the pace of evolution would be very slow, because a population would not attain a total population size that would make it probable that rare beneficial mutations would occur and become predominant in a population.

Even though it was mentioned that genetic evidence suggests that the present world human population are all descendants of a small band of our species that migrated out of Africa about 60,000 years ago, little attention in this paper was given to the evolution of the world population over time. In this connection, the book by McEvedy and Jones [21] would be of interest in which estimates of population size are given from about 400 years before the common era (BCE) to 1975. Another book of interest is that of McKeown [22], which traces the evolution of the European population, and in particular the population of England and Wales, where the rise modern population did not start to occur until about 1750. As pointed out in this book, the evolution of modern medicine played a significant role recent rise of population. The book by Hostetler [23] is also of interest, because it provides a documentation of the Hutterite society over a little more than a century, since they migrated to South Dakota, USA, in the 19th century and have expanded their population to various communities in Montana and in neighboring prairie provinces of Canada. The Hutterites are of particular interest, because, as a population with records spanning over a century, they experienced one of the highest levels of fertility in recorded human history.

As this paper is being written, the estimated size of the world human population is at about 7 billion, and as this population grows, a question that naturally arises is how many people can the earth support? In a definitive book with this title, Cohen [24] addresses this question in depth. Among the many issues addressed in this book are questions regarding the carrying capacity of the earth. It is very difficult to provide usable answers to this question, but it may turn out that the approach used to formulate the notion of a carrying capacity suggested in this paper for an age structured population may be useful in obtaining experimental answers to this question. But it should also be mentioned that other approaches yet to be invented may also be useful in an experimental approach to estimate the carrying capacity of the earth.

It should also be mentioned that climate change has also played an important role in human evolution, and moreover, catastrophic events during the past have led to an abrupt climate change. An example of such a change was the eruption of the mega Toba volcano on the island of Sumatra about 70,000 to 75,000 years ago. According to Rampino and Ambrose [25] during this mega eruption vast quantities of materials emanating from the earth’s crust were thrust into the earth’s atmosphere resulting in a “nuclear winter” that cooled the earth to the extent that there were frosts in the tropics for five to ten years that resulted in population crash of animals and plants, including populations of humans. According to this view, all human on the earth today are descendants of a relatively small population that survived the Toba eruption. A common term used by evolutionists to describe small population size during evolution of a species is “bottleneck.”

The short list of references cited after that are informative but not current. It is, therefore, of interest to search for more recent publications on the evolution of world population. When this phrase was typed into an internet search engine, a total of about 29,100,000 sites were listed. Within this large listing of sites, a large number of technical papers were also listed that were beyond the scope of this paper. But, there were two significant papers that were within the scope of this paper. One on these papers was that of Hawks et al. [26] in which evidence for a population bottleneck was investigated during the Pleistocene period of human evolution that included the Toba eruption. Another paper was that of Stearns et al. [27] in which the measurement of selection was investigated, using data from contemporary populations. In comparison with an average human lifespan, evolution of a population occurs over very lager periods of time expressed in terms of generations. Therefore, when a contemporary set of data is used in attempts to measure selection, it is necessary to investigate the heritability of a trait measured from generation to generation, and as this paper is read it is recommended that a reader be mindful of the concept of heritability. For if an expression genotype or phenotype is not governed by some region (or regions) of a genome, the trait cannot be passed on from generation to generation and thus cannot be subject to natural selection.

In the evolving literature on evolutionary and population genetics, two themes have been received considerable attention, namely, genome wide searches for signatures of selection and methods for simulating models of genomes. In genome wide searches, a problem that often arises is that of developing statistical criteria to differentiate among regions of a genome, such that selection has acted in the recent evolutionary past from regions that presumedly have evolved under neutral evolution. Two papers describing methods for detecting signatures of positive selection are those of Sabeti et al. [28] and Grossman et al. [29]. Additional references on this subject may be found in the cited literature in Mode and Gallop [15]. As yet methods have not been developed to model the evolution region of a genome with respect to mutation and selection within the framework of the age structured process under consideration, but if a reader is interested in a brief review of published papers dealing with simulating the evolution of a genomic region, the papers cited in chapter 14 of Mode and Sleeman [9] may be consulted. Among these papers there are Chadeau-Hyam et al. [30], Hoggart et al. [31], and Hoggart et al. [32]. Preliminary algorithms for simulation various types of mutations, such as nucleotide substitutions, deletions, insertions, and repeats are also included in chapter 14 of Mode and Sleeman.

References

B. Charlesworth, Evolution in Age-Structured Populations, vol. 1 of Cambridge Studies in Mathematical Biology, Cambridge University Press, Cambridge, Mass, USA, 1980.

T. E. Harris, The Theory of Branching Processes, Die Grundlehren der Mathematischen Wissenschaften, Bd. 119, Springer, Berlin, Germany, 1963.

K. B. Athreya and P. E. Ney, Branching Processes, Springer, New York, NY, USA, 1972.

P. Jagers, Branching Processes with Biological Applications, John Wiley & Sons, London, UK, 1975.

C. J. Mode, Multitype Branching Processes: Theory and Applications, American Elsevier, New York, NY, USA, 1971.

S. Asmussen and H. Hering, Branching Processes, vol. 3 of Progress in Probability and Statistics, Birkhäuser, Boston, Mass, USA, 1983.

M. Kimmel and D. E. Axelrod, Branching Processes in Biology, vol. 19 of Interdisciplinary Applied Mathematics, Springer, New York, NY, USA, 2002.

P. Haccou, P. Jagers, and V. A. Vatutin, Branching Processes: Variation, Growth, and Extinction of Populations, Cambridge Studies in Adaptive Dynamics, Cambridge University Press, Cambridge, UK, 2007.

C. J. Mode and C. K. Sleeman, Stochastic Processes in Genetics and Evolution: Computer Experiments in the Quantification of Mutation and Selection, World Scientific Publishing, Hackensack, NJ, USA, 2012.

C. J. Mode, Stochastic Processes in Demography and Their Computer Implementation, vol. 14 of Biomathematics, Springer, Berlin, Germany, 1985.

C. J. Mode and C. K. Sleeman, Stochastic Processes in Epidemiology, HIV/AIDS, Other Infectious Diseases and Computers, World Scientific, Singapore, 2000.

C. J. Mode, T. Raj, and C. K. Sleeman, “Monte Carlo implementations of two sex density dependent branching processes and their applications in evolutionary genetics,” in Applications of Monte Carlo Methods in Biology, Medicine and Other Fields of Science, C. J. Mode, Ed., pp. 273–296, INTECH, 2011.

C. J. Mode, T. Raj, and C. K. Sleeman, “Simulating the emergence and survival of mutations using a self regulating multitype branching processes,” Journal of Probability and Statistics, vol. 2011, Article ID 867493, 39 pages, 2011.

J. Hawks, The Rise of Humans: Great Scientific Debates, The Teaching Company, Chantilly, Va, USA, 2011, http://www.thegreatcourses.com/.

C. J. Mode and R. J. Gallop, “A review on Monte Carlo simulation methods as they apply to mutation and selection as formulated in Wright-Fisher models of evolutionary genetics,” Mathematical Biosciences, vol. 211, no. 2, pp. 205–225, 2008.

L. Deng and J. Lin, “Random number generation of a new century,” The American Statistician, vol. 54, no. 2, pp. 145–150, 2000.

P. L'Ecuyer and R. Touzin, “On the Deng-Lin random number generators and related methods,” Statistics and Computing, vol. 14, no. 1, pp. 5–9, 2004.

L. Deng, “Efficient and portable multiple recursive generators of large order,” ACM Transactions on Modeling and Computer Simulation, vol. 15, no. 1, pp. 1–13, 2005.

M. Alderson, International Mortality Statistics, Facts on File, Inc., New York, NY, USA, 1981.

D. L. Hartl and A. G. Clark, Principles of Population Genetics, Sinauer Associates, INC., Sunderland, Mass, USA, 1989.

C. McEvedy and R. Jones, Atlas of World Population History, Facts on File, New York, NY, USA, 1979.

T. McKeown, The Modern Rise of Population, Academic Press, New York, NY, USA, 1976.

J. Hostetler, Hutterite Society, The Johns Hopkins University Press, Baltimore, Md, USA, 1974.

J. E. Cohen, How Many People Can the Earth Support?W.W. Norton & Company, New York, NY, USA, 1995.

M. Rampino and S. Ambrose, “Volcanic winter in the Garden of Eden: the Toba supereruption and the late Pleistocene human population crash,” in Symposium: Catastrophism, Natural Disasters and Cultural Change, World Archaeological Congress 4, University of Cape Town, 1999.

J. Hawks, K. Hunley, S. Lee, and M. Wolpoff, “Population Bottlenecks and Pleistocene evolution,” Molecular Biology and Evolution, vol. 17, no. 1, pp. 2–22, 2000.

S. C. Stearns, S. G. Byars, D. R. Govindaraju, and D. Ewbank, “Measuring selection in contemporary human populations,” Nature Reviews Genetics, vol. 11, pp. 611–622, 2010.

P. Sabeti, et al., “Positive selection in the human linage,” Science, vol. 312, pp. 1614–1620, 2006.

S. Grossman, et al., “A composite of multiple signals distinguishes causal variants in regions of positive selection,” Science, vol. 327, pp. 883–886, 2010.

M. Chadeau-Hyam, et al., “Fregene: simulation of realistic sequencelevel-data in populations and ascertained samples,” BMC Bioinformatics, vol. 9, article 364, 2008.

C. Hoggart, et al., “Sequence level population simulations over large genomic regions,” Genetics, vol. 177, pp. 1725–1731, 2007.

C. Hoggart, et al., “Genome-wide significance for dense SNP and resequencing data,” Genetic Epidemiology, vol. 32, pp. 179–185, 2008.