Review Article | Open Access

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

**Academic Editor:**Peter Olofsson

#### 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.

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 .

As can be seen from this figure, significant numbers of individuals of mutant genotype