Markov chain approach to the distribution of ancestors in species of biparental reproduction
MMarkov chain approach to the distribution of ancestors in species of biparental reproduction
M. Caruso
Departamento de F´ısica Te´orica y del Cosmos, Universidad de Granada,Campus de Fuentenueva, Granada (18071), Espa˜na
C. Jarne
Departamento de F´ısica, Facultad de Ciencias Exactas, IFLP-CONICET,Universidad Nacional de La Plata, La Plata (1900), C.C.67, Argentina
We studied how to obtain a distribution for the number of ancestors in species of sexual repro-duction. Present models concentrate on the estimation of distributions repetitions of ancestors ingenealogical trees. It has been shown that is not possible to reconstruct the genealogical history ofeach species along all its generations by means of a geometric progression. This analysis demon-strates that it is possible to rebuild the tree of progenitors by modeling the problem with a Markovchain. For each generation, the maximum number of possible ancestors is different. This bringshuge problems for the resolution. We found a solution through a dilation of the sample space, al-though the distribution defined there takes smaller values respect to the initial problem. In order tocorrect the distribution for each generation, we introduced the invariance under gauge (local) groupof dilations. These ideas can used to study the interaction of several processes and provide a newapproach on the problem of the common ancestor. In the same direction, this model also providesome elements that can be used to improve models of animal reproduction.
INTRODUCTION
Up till now, previous attempts aiming to calculate thenumber of ancestors in species of sexual reproductionhave not been totally successful. Present models con-centrate on the estimation of distributions of ancestorsrepetitions in genealogical trees [1–3]. It has been shownthat is not possible to reconstruct the genealogical his-tory of each species along all its generations by meansof a geometric progression [4]. The reason for that isthe geometric progression is determined by a sequenceof independent events. We postulate that blood relation-ship is a kind of interaction that connects the events. Itis possible to re-build the tree of progenitors by mod-elling the problem with a Markov chain. If we consider arandom variable which represents the number of ances-tors present in a given generation, the size of the samplespace depends of each generation. This brings seriouscomplications on the solution of the problem, not only ofmathematical nature. We propose submerge the originalsample space into a larger one. This dilution modifies the probability distribution. We show the need to imple-ment a covariant derivative, due to a gauge transforma-tion, which leaves the evolution equation invariant andcorrect the probability distribution.The main goal of present work is describe the distri-bution of ancestors for species with sexual reproduction,but also show the novel method used here to solve otherstochastic problems.There are two important assumptions about the bi-ology of the considered system. The first one is thatthe species described here has not specific behavior ofsexual partner selection (random mating reproduction)[5, 6]. Many species or population groups exhibit thiskind of reproduction. This is the most simple case toperform the calculation. The second assumption is aboutthe population size. The distribution of ancestors fora given generation is contained in a population largeenough to not force the selection of sexual partners bloodrelated. The partners could be blood related or not,randomly. Current model presents a random matingin non-overlapping generations with negligible mutation a r X i v : . [ phy s i c s . b i o - ph ] O c t and selection. These two assumptions are common to de-velop population genetic models, in particular, these arepresent in the Hardy-Weinberg principle [7, 8].In this work we show a way to calculate a probabilitydistribution to get a certain number of ancestors for eachgeneration. We have obtained its first two cumulants:the expected value of the number and its dispersion.Specific conditions about small size populations, orspecific sexual behavior can be considered later as mod-ifications of the general case described here.Present work will be useful in order to understand theorigin of species extrapolating the individual genealogyfor all members at the beginning of the species. It is pos-sible to go one step further, to establish how populationscan be affected by certain conditions, such as isolation ormigration of individuals, by studying population groupswith different genetic pool [9]. These ideas can be used toperform more realistic models in animal populations andalso, improve estimations about extinction processes. A MARKOVIAN APPROACH TO THEANCESTORS PROBLEM
To calculate the number of ancestors of an individ-ual it is necessary to use a statistical approach. If wesimply accept that 2 t +1 allow us to calculate the num-ber of ancestors in the t − generation, where t = 0 is thegeneration of progenitors of the first order (or parents forshort) and so forth, we arrive at an absurdity. As we turnto past generations, the probability that some ancestorshave been relatives is significantly larger [1–3]. This im-plies a restriction on the number of ancestors with respectto 2 t +1 . This last quantity corresponds to the maximumpossible number of ancestors in each t − generation.There are several examples showing different ways inwhich the number of ancestors is reduced with respectto the maximum number in each generation. As statedin [4], the reduction of the number of ancestors, com-pared to 2 t +1 , is caused by blood relationship. Figure 1shows, as an example, only the three first generations ofgenealogical tree with two different ways to constrain thenumber of ancestors. There is a way to weight the bloodrelationship using a statistical approach that includes allpossible kinds of relationship in each generation. In this approach the only constraint in the number of ancestorsis caused by random blood relationship between individ-uals of the same generation. We considered a populationof ancestors whose maximum size in each t − generationis given by the geometric progression 2 t +1 . FIG. 1. (Color online) Examples of three kinds of genealogicaltrees (only the first few g − generations) a : No restrictions byblood relationship. b and c : Two kind of restriction in thirdgeneration, ancestors sharing one (b) or two parents (c) . Therestriction by blood relationship increases according to thedegree of endogamy. We did not consider any restriction for the numberof ancestors generated by issues related to culture, in thehuman case, ethological in the animal case, or isolation ofpopulations, etc. If we want to study the distributions ofancestors of individuals from populations where there areless or equal individuals than 2 t +1 for the t − generation,there is an additional restriction on the number of ances-tors. Blood relationship interconnects the events in theoriginal process that leads to 2 t +1 , which was generatedby independent events and no relation between ancestorsof each generation.Derrida’s model [1, 2] is based on numerical simulationunder the same assumptions (closed population evolvingunder sexual reproduction with non-overlapping genera-tions). The population size is fixed for all generationsand equally divided into two groups, representing malesand females. At every generation, they form randomheterosexual pairs and assign them a certain number ofdescendants according to a Poisson distribution. Thisis done by choosing for each male or female one pair ofparents at random in the previous generation.In our work the population size is not fixed, but alwaysbigger than 2 t +1 for each t − generation.We defined two random variables y ( t ) and x ( t ) whichrepresents the number of individuals who are inside andoutside to the set of ancestors, with respect to the maxi-mum possible number of ancestors in each t − generation.For this definition we have x ( t ) + y ( t ) = 2 t +1 . (1)We considered each generation as a link of the chainwhich form a first order Markov process [10, 11]. Thisprocess is constructed on a given set of individuals or-dered by generations. We take the current generationand we count its parents. Then we take all selected in-dividuals and remake the previous question and so forth.There exists a generation in which the question or pre-vious classification makes no more sense, in which casethe process ends after a given generation. This kind ofprocess is widely used to describe the evolution of traitsthat adopt a finite number of states [12].From equation (1) we have y ( t ) = 2 t +1 − x ( t ) , (2)if y ( t ) describes a Markov process that implies x ( t ) de-scribes another Markov process.We do not distinguish the different kinds of blood rela-tionship between the ancestors of a particular generationsuch as brothers or cousins, and so forth. We simply con-sider them as indistinguishable and we just count howmany there are. For the purpose of the calculations weconsider t as a continuous variable. Finally we associatea discrete-time Markov process to the continuous-timeMarkov process { x ( t ) : t ≥ } called a skeleton process [13] defined as { x ( g ) : g ≥ } , where g is the generationnumber. The time evolution of this process is determined bythe knowledge of the probability distribution in each t − generation, denoted by p n ( t ) = P [ x ( t ) = n ] (3)for all ( n, t ) ∈ S t × R , where S t is the sample space of x ( t )which corresponds to the interval [0 , n t ], n t = (cid:98) t +1 (cid:99) − (cid:98) z (cid:99) is the integer part of a real number z .An equivalent way to describe the process it is throughan initial value p n (0) and the conditional probabilitygiven by P nm ( t, s ) = P [ x ( t ) = n | x ( s ) = m ], whichrepresents the transition matrix elements of the states( m, s ) (cid:55)−→ ( n, t ).For each generation the events are mutually exclusive.Consequently at the time t + (cid:15) the probability of find n re-strictions is given by to the transition from m restrictionsat the time t , in this way p n ( t + (cid:15) ) = (cid:88) m ∈ S t P nm ( t + (cid:15), t ) p m ( t ) . (4)After some elementary operations (see appendix section A1 ) we get d t p n ( t ) = (cid:88) m ∈ S t Q nm ( t ) p m ( t ) (5)where d t denotes the total time derivative ddt , Q nm ( t ) = ∂ t P nm ( t, s ) | s = t is called the infinitesimal generator and δ nm is the Kroneker delta.We define ϕϕϕ ( t ) as an | S t | -tuple of the probability dis-tribution as ϕϕϕ ( t ) = ( p ( t ) , p ( t ) , · · · , p n t ( t ) ) (cid:124) , where | S t | denotes the cardinal number of S t and (cid:124) represents thetransposition.The evolution equation for the process can be ex-pressed in a matrix form as d t ϕϕϕ ( t ) = QQQ ( t ) ϕϕϕ ( t ) (6)We denote the expectation number of ancestors by α ( t ) = (cid:104) y ( t ) (cid:105) and from the equation (2) α ( t ) = 2 t +1 − (cid:104) x ( t ) (cid:105) (7)where (cid:104) x k ( t ) (cid:105) is the expectation value of x ( t ) raised tothe positive integer power k (or k -moment for short)of the distribution p n ( t ) and by definition is (cid:104) x k ( t ) (cid:105) = (cid:80) n n k p n ( t ). The quantity (cid:104) x ( t ) (cid:105) represents a constraintcaused by blood relationship, which affects the expecta-tion number of ancestors in each generation. DILUTION OF SAMPLE SPACE VIA GAUGEGROUP OF DILATIONS
The sample space of x ( t ) is different for each t − generation, thus there is enormous difficulty to solvethe equation (6). We considered a dilution of S t withina larger set S ⊇ S t , for all t , consisting of replacing theendpoint n t by a huge number N . This dilution can beviewed as a dilation represented with the substitutionrule n t (cid:55)−→ N , such that S = [0 , N ]. On the other handwe know that there exist a certain T − generation that canbe considered as the end of the process. The existence ofa limit generation, T , allows us to choose N = n T . Con-sequently we can solve the problem in this dilated samplespace and then recover the lost endpoint caused by thedilation through a suitable transformation. The price topay for it is the need of renormalization of the distribu-tion defined on S to compensate the dilution effect. Therenormalization takes place by a linear transformationwhich modifies the norm of the distribution for each gen-eration. This local transformation (i.e. depends of each t ), is structured as a gauge group , specically the groupof local dilations . Essentially, the distribution defined on S t is equivalent to the renormalized distribution which isdefined on the dilated sample space S .In summary, we can interpret that the process on S t is the result of a process on this larger set S which in-teracts with another process on the complement of S t ,denoted by S − S t . This interaction is represented bythe renormalization of the distribution defined on S , inan effective theory context. As long as the process on S becomes much simpler, the description on S − S t will bemore complex. This is the basis for the dilation transfor-mation, which is discussed in the appendix sections A5 and A6 .We considered a version in which the sample space S t is dilated to the set of natural numbers N , including the0 element. Then we have only one boundary conditionfor the state n = 0. This allows us to focus on timehomogeneous processes , i.e. the infinitesimal generatoris independent of t . Another consideration is the spatialhomogeneity , i.e. the case where the infinitesimal gener-ator does not depend on the state of the random variable X ( t ). The Markov process in this larger sample space N re-quires to consider two new random variables { X, Y } de-fined on N and related in a similar way to the old randomvariables { x, y } from (2). The associated probability dis-tribution is denoted by P n ( t ) = P (cid:2) X ( t ) = n (cid:3) and defines φφφ ( t ) = ( P ( t ) , P ( t ) , · · · ) (cid:124) which satisfies the equation d t φφφ ( t ) = Q φφφ ( t ) . (8)Knowing the initial conditions φφφ (0) = (1 , , · · · ) (cid:124) andthe infinitesimal generator Q we can write the formalsolution of (8) as φφφ ( t ) = exp( t Q ) φφφ (0) . (9)In order to establish the matrix Q , we study the timeevolution t (cid:55)−→ t + (cid:15) for small value of (cid:15) . Therefore, onlytransitions to the nearest states are allowed, because theinfinitesimal time evolution only has a finite variety oftransition states. For n (cid:54) = 0 these transitions are n (cid:55)−→{ n − , n, n + 1 } and n (cid:55)−→ { n, n + 1 } , for n = 0.Considering this brief discussion, the dynamics de-scribed by the equation (8) and the imposed conditionsrepresents a time homogeneous birth-death process . Anaive way to picture the process in the context of queue-ing theory [14], is through one queue representing all an-cestors waiting to be classified if they are blood relatedor not by one server.In the appendix section A4 , we show how to choosea numerical matrix Q . Finally the evolution equationstakes the form d t P n ( t ) = P n +1 ( t ) − P n ( t ) + P n − ( t ) , (10) d t P ( t ) = P ( t ) − P ( t ) , together with the initial condition which is P n (0) = δ n ,we obtain the explicit solution [14] P n ( t ) = e − t [ I n (2 t ) + I n +1 (2 t )] (11)where I n ( x ) is the modified Bessel function [15]. A briefdescription to obtain the solution (11) is also present in[14]. There is a construction of the generatrix function g ( t, z ) = (cid:80) n ∈ N P n ( t ) z n , see (A.22), and from (8) derivean equation for g ( t, z ).The equation (10) is the generic expression for allMarkov processes on denumerable sample spaces andcontinuous time with a particular values of Q . THE GAUGED DISTRIBUTION OF ANCESTORS
As we have previously argued, before using this dis-tribution to calculate the moments, it is necessary toperform a renormalization process. The reason is thatthe solution given by (11) is normalized over N . We per-form a gauge transformation [16], denoted by g t , whichis applied to the probability distributions as g t : P n ( t ) −→ λ ( t ) P n ( t ) . (12)The transformation (12) leaves the evolution equation(8) invariant and allows both distributions to describea Markov process. We denote p n ( t ) = λ ( t ) P n ( t ) thegauge transformed distribution of P n ( t ). The action ofthe group g t applied to the distribution P n ( t ) leads toa distribution p n ( t ) defined over S t . This idea can beunderstood in the context of conditional probabilities,with which we can obtain a projection of the distributionon N into S t , keeping the correct normalization.In oder words, the transformation g t leads to a newrandom variable X , which is the gauge transformed of X .To preserve the invariance of (8) under g t , we introducea covariant derivative D t = d t − ω ( t ) (13)where ω ( t ) = d t λ ( t )[ λ ( t )] − . See the appendix for a moreextensive explanation.The expectation value of X raised to a positive integerpower k is (cid:104) X k ( t ) (cid:105) = (cid:80) n n k p n ( t ). This allows us to writea general relation between (cid:104) X k ( t ) (cid:105) and (cid:104) X k ( t ) (cid:105)(cid:104) X k ( t ) (cid:105) = λ ( t ) (cid:104) X k ( t ) (cid:105) . (14)Rescaling the process described by X ( t ) and using thesolution (11) we calculated the first two cumulants (cid:104) X ( t ) (cid:105) = λ ( t ) (cid:104) X ( t ) (cid:105) , (15) (cid:104) [ X ( t ) − (cid:104) X ( t ) (cid:105) ] (cid:105) = λ ( t )[2 t − (cid:104) X ( t ) (cid:105) − (cid:104) X ( t ) (cid:105) ] . where (cid:104) X ( t ) (cid:105) = e − t (cid:2) t I (2 t ) + (cid:0) t + (cid:1) I (2 t ) (cid:3) − . (16)From equation (2) the variance of x is equal to thevariance of y . The same argument is valid for X and Y . We define the standard deviation of Y ( t ) = 2 t +1 − X ( t ),denoted by σ ( t ), as the square root of the second equationof (15), which quantifies the statistical fluctuation.As we consider a constant function ω ( t ), then λ ( t ) = 2 a t + b . (17)We have obtained a family of functions for the expec-tation number of ancestors α ( t ) = 2 t +1 − λ ( t ) (cid:104) X ( t ) (cid:105) (18)parametrized by the real numbers a and b of (17).If the expected value satisfies α ( t ) = α and α ( t ) = α , for two generations t and t such that t (cid:54) = 0 (cid:54) = t ,the parameters a and b can be obtained by a = 1 t − t log (cid:20) t +1 − α t +1 − α (cid:104) X ( t ) (cid:105)(cid:104) X ( t ) (cid:105) (cid:21) (19) b = 1 t − t (cid:40) t log (cid:20) t +1 − α (cid:104) X ( t ) (cid:105) (cid:21) − t log (cid:20) t +1 − α (cid:104) X ( t ) (cid:105) (cid:21) (cid:41) where naturally α i ≤ t i +1 , for i = 1 ,
2, to ensure gooddefinition of a and b .The gauge transformation modulates the amplitude of (cid:104) X ( t ) (cid:105) . This allows us to define the notion of horizon-tal and vertical range of α ( t ). One important point ofthe curve α ( t ) is the maximum generation range , this isa nonzero generation T in which α becomes equal to 2.Another interesting point is the maximum of α ( t ), whichdetermines the intensity of the process . Without loss ofgenerality we can choose t = T , in which case α ( t ) = 2,and α ( t ) = sup { α ( t ) : t ∈ [0 , T ] } . For any pair of differ-ent points, ( t , α ) and ( t , α ), considered relevant, weselect one and only one curve of the family, parameterizedby a and b given by (19). The gauge transformation g t ,through the a and b parameters, controls both horizontalrange and vertical range of the process. This T maybenot be a realistic value, but fix a maximum number ofgenerations of a particular species may have.For illustrative purposes, in figure 2, we have selected 4curves to the expectation number of ancestors α ( t ) givenby (18) and parametrized by different values of { a , b } .We include a geometric progression 2 t +1 , which corre-sponds to the maximum possible number of ancestors ineach t − generation. FIG. 2. (Color online) The 4 values of { a , b } are ob-tained from (19) and parametrized from ( t , α , t , α ) =( t , ξ t +1 , T, ≤ ξ ≤
1. In t the curve reachesthe fraction ξ of the total number of possible ancestorsfor that generation, while t defines the maximum genera-tion range denoted by T . The geometric progression 2 t +1 is in dash black line. The { a,,,b,,,c,,,d } ( { blue ,,, green ,,, orange ,,, red } ) lines can be obtained respectively from ( t , α , t , α ) ∈ (cid:8) (3 , . × , , ,,, (3 , . × , , ,,, (3 , . × , , ,,, (3 , . × , , (cid:9) . Figure 3 shows three realizations of the number of an-cestors in terms of the expectation value α ( t ) and a mea-sure of the dispersion given by σ ( t ), for a particular valuesof a and b .This model may be employed in order to recognizea possible threshold to identify high endogamic popu-lations as well as its possible causes. Using the genealog-ical tree, the model can be used to indicate which livingspecies may be near to the extinction. FINAL COMMENTS AND POSSIBLE MODELEXTENSIONS
The model explained above allows to calculate the ex-pectation number of ancestors in each generation, consid-ering the possibility of blood relationship between indi-vidual of the same generation and a population of ances-tors which maximum size is 2 t +1 . But there are two pos-sible generalizations. The model can be extended to takeinto account relationships between individuals of otheradjacent generations using a similar idea, simply consid-ering higher order Markov chains. By introducing the FIG. 3. (Color online) A band of curves, defined by the set B = { ˆ α ( t ) : ˆ α ( t ) ∈ [ α ( t ) − σ ( t ) , α ( t ) + σ ( t )] } , contains the ex-pectation value α ( t ) given by the equation (18) and its statis-tical fluctuation at 1- σ . The { a,,,b,,,c } ( { blue ,,, green ,,, red } ) linescorresponds to the realization of { α ( t ) − σ ( t ) , α ( t ) , α ( t )+ σ ( t ) } ,with the values of ( t , α , t , α ) = (6 , . × , , corresponding terms in the infinitesimal generator Q , forexample an absorbing barrier [17], the same model canbe used to calculate the expectation number of ancestorsin a specific population with additional restrictions suchas isolation, immigration, specific reproductive behavior,or cultural restrictions for the human case.In future applications we could generalize the modelthrough a new constraint to fix the maximum number ofindividuals at certain generation. This proposal implies ageneralization of this work in which the maximum num-ber of ancestors will be given by a piecewise function γ ( t )instead of 2 t +1 . This leads to slightly modify the pro-cess defined at the beginning in (2) as y ( t ) = γ ( t ) − x ( t )and the endpoint of the sample space for x ( t ) be comes n t = γ ( t ) −
2. This generalization includes a time inho-mogeneity in the infinitesimal generator Q and preservesan appropriate renormalization.We can include these possible extensions using the pro-cess { X ( t ) } and study a most general gauge g t transformgiven by a linear transformation g t : { P m ( t ) } m ∈ N −→ p n ( t ) = (cid:88) m ∈ N λ nm ( t ) P m ( t ) (20)where λλλ = { λ nm } is a non singular matrix. We canexpress the last expression in a matrix form g t : φφφ −→ φφφ (cid:48) = λφλφλφ (21)where φφφ (cid:48) ( t ) = ( p ( t ) , p ( t ) , · · · ) (cid:124) .In order to preserve the invariance of (8) under thisgeneralization of g t , we introduce the corresponding co-variant derivative D t = d t − ωωω ( t ) (22)where ωωω ( t ) = d t λλλ ( t )[ λλλ ( t )] − .The evolution equation for φφφ (cid:48) ( t ) is also invariant underthe local dilation group. D t φφφ (cid:48) ( t ) = Q (cid:48) ( t ) φφφ (cid:48) ( t ) (23)and the gauged infinitesimal generator is now Q (cid:48) ( t ) = λλλ ( t ) Q [ λλλ ( t )] − (24)which corresponds to a similarity transformation of Q .This model itself can be applied to describe other bi-ological or physical systems with similar dynamics. Sta-tistical models of biparental reproduction have alreadybeen compared with physical systems before, such as,spin-glass systems [18]. In this regard the evolutionarygraph theory is an approach to study how topology affectsthe evolution of a population [19].Other analogous processes to the biparental reproduc-tion in physics are described with similar statistical ormarkovian models [10]. In high energy physics the pro-duction of a cascade by a cosmic ray is described by theHeitler model [20]. Although this model is different fromthe one presented here we could compare the number ofancestors with the number of particles in each genera-tion and reinterpret this results in terms of these kindsof phenomena.It is possible to estimate the maximum generationrange, T , searching in the fossil record the first time thata particular species appears and use its reproductive rate.In this way we are classifying each species not in termsof life time on earth (time units), but according to thenotions of generational patterns.The interaction of various of these processes can becombined with universal common ancestor’s models [21]to understand the development of a certain species. The ideas in the current model can be used in biology, pop-ulation ecology and genetics. An important achievementof the model is that based on the previous knowledge ofthe life time of a certain species, we can calculate thenumber of ancestors in each generation of this species.More fundamental uses of this ideas can be found inmathematics, in particular in theory of stochastic pro-cesses and physics area connected with the theory ofstochastic processes. Future research through a La-grangian description may find novel applications of thepresent proposal. In this case we will consider the prob-abilities { p n ( t ) } as the set of generalized coordinates. ACKNOWLEDGMENTS
We thank our respective PhD advisors: Fernando Cor-net and Mar´ıa Teresa Dova, Hern´an Wahlberg. Also wethank to Carlos Garc´ıa Canal and Huner Fanchiotti fortheir helpful criticism. We are indebted to Federico Ag-nolin and Ben Page for reading the manuscript and pro-viding advice. Recall also our anonymous readers andreviewers for their contribution to this work. And fi-nally a special mention to Micaela Moretton, Mar´ıa ClaraCaruso and Gabriel Lio for local support.
APPENDIXA1
ON THE EVOLUTION EQUATION IN S We defined the random variable y ( t ) associated withthe number of ancestors as y ( t ) = 2 t +1 − x ( t ) (A.1)where S t is the sample space of x ( t ) and t is a continuousvariable. Then, according to a discretization process, thedistribution of ancestors is obtained and the variable t will be the number of generations.We denoted by α ( t ) the expectation number of y ( t ) as α ( t ) = 2 t +1 − (cid:104) x ( t ) (cid:105) (A.2)where (cid:104) x ( t ) (cid:105) constrains the number of ancestors underthe blood relationship hypothesis. We will focus on theprocess of this random variable x ( t ) and our probabilitydistribution denoted by p n ( t ) = P [ x ( t ) = n ].In the problem of counting ancestors the sample space S t is different for each t -generation. For this reason weintroduced a fictitious process on a larger sample space S than the original S t . Therefore we have solved theproblem in S and renormalized the distribution to takeinto account of the interaction with the lost boundary of S t . In other words, we have considered a dilution of S t into S . In particular we used the set of natural numbers S = N , including the 0 element.We considered the study of the evolution over a genericdenumerable sample space S = [0 , N ] ⊂ N with the ran-dom variable X ( t ) and probability distribution P n ( t ) = P [ X ( t ) = n ]. This evolution is governed by the condi-tional probability given by P nm ( t, s ) = P [ X ( t ) = n | X ( s ) = m ] (A.3)The matrix P ( t, s ) satisfies the Chapman-Kolmogorov equation [11, 17] P ( t, s ) = P ( t, u ) P ( u, s ) (A.4)for 0 ≤ s ≤ u ≤ t . Also the sum of the elements of eachcolumn is (cid:88) n ∈ S P nm ( t, s ) = 1 . (A.5)For the general case we develop a power series of thematrix P ( t + (cid:15), s ), for a fixed value of s , we have P nm ( t + (cid:15), s ) = P nm ( t, s ) + (cid:15) ∂ t P nm ( t, s ) + . . . (A.6)where ∂ t is a simplified notation of partial time derivative ∂∂t In order to obtain the equation (4) we study the timeevolution t (cid:55)−→ t + (cid:15) , for small value of (cid:15) . We need toknow P nm ( t + (cid:15), t ) then (A.6) becomes P nm ( t + (cid:15), t ) = δ nm + (cid:15) ∂ t P nm ( t, s ) | s = t + . . . (A.7)we recognize the second term of (A.8) as the infinitesimalgenerator Q nm ( t )Q nm ( t ) = lim (cid:15) → P nm ( t + (cid:15), t ) − δ nm (cid:15) (A.8) We assumed that t is continuous. This allows us toevaluate the process at any t between two generations,but it also reduces the number of possible states in aninfinitesimal evolution. Therefore only transitions to thenearest states are allowed. For the ancestry problem,the infinitesimal time evolution has a finite number oftransition states. These transitions are denoted by n (cid:55)−→ n (cid:48) , where n (cid:48) ∈ T n , i.e. n (cid:48) depends on the initial state n .We write explicitly T = { , } , T N = { N − , N } andfor n (cid:54) = 0 , N : T n = { n − , n, n + 1 } . Furthermore, if | n (cid:48) − n | > P ( t + (cid:15), t ) is normalizedfor all t and small (cid:15) as (cid:88) n (cid:48) ∈ T n P n (cid:48) n ( t + (cid:15), t ) = 1 (A.9)note that n (cid:48) runs over T n , depending on whether n isequal to 0, N or any other value of S − { , N } .We express P n (cid:48) n ( t + (cid:15), t ) for these three cases from(A.8) P n − n ( t + (cid:15), t ) = µ n ( t ) (cid:15) + O t ( (cid:15) ) , P n +1 n ( t + (cid:15), t ) = ν n ( t ) (cid:15) + O t ( (cid:15) ) , (A.10) P n (cid:48) n ( t + (cid:15), t ) = 0 , | n (cid:48) − n | > , P nn ( t + (cid:15), t ) = 1 − [ ν n ( t ) + µ n ( t )] (cid:15) + O t ( (cid:15) ) , where O t ( x ) represents a type of function that goes tozero with x faster than x , for a given t , that islim x → O t ( x ) x = 0 . (A.11)The fourth equation of (A.10) is obtained through thefirst three of them. In the general case we proceed asfollows from (A.5) P n n ( t, s ) = 1 − (cid:88) n (cid:48) ∈ S −{ n } P n (cid:48) n ( t, s ) (A.12)then from an infinitesimal time evolution s ≡ t (cid:55)−→ t + (cid:15) and (A.9) P nn ( t + (cid:15), t ) = 1 − (cid:88) n (cid:48) ∈ T n −{ n } P n (cid:48) n ( t + (cid:15), t ) P nn ( t + (cid:15), t ) = 1 − P n +1 n ( t + (cid:15), t ) − P n − n ( t + (cid:15), t )which is the fourth equation of (A.10).Replacing (A.10) in (4) and written for n (cid:54) = 0, we have P n ( t + (cid:15) ) = [ ν n − ( t ) (cid:15) + O t ( (cid:15) )] P n − ( t )+ [ µ n +1 ( t ) (cid:15) + O t ( (cid:15) )] P n +1 ( t )+ [1 − ν n ( t ) (cid:15) − µ n ( t ) (cid:15) + O t ( (cid:15) )] P n ( t )then P n ( t + (cid:15) ) − P n ( t ) (cid:15) = (cid:20) ν n − ( t ) + O t ( (cid:15) ) (cid:15) (cid:21) P n − ( t )+ (cid:20) µ n +1 ( t ) + O t ( (cid:15) ) (cid:15) (cid:21) P n +1 ( t ) − (cid:20) ν n ( t ) + µ n ( t ) + O t ( (cid:15) ) (cid:15) (cid:21) P n ( t )taking the limit (cid:15) → d t P n ( t ) = ν n − ( t ) P n − ( t ) + µ n +1 ( t ) P n +1 ( t )(A.13) − [ ν n ( t ) + µ n ( t )] P n ( t )The stochastic process described by equation (A.13)corresponds to the general class of stochastic dynamicscalled birth and death process , which includes the queue-ing process [11, 14].The functions µ n ( t ) and ν n ( t ) are part of the infinites-imal generator Q ( t ). From the first equation of (A.10)we have µ n ( t ) = lim (cid:15) → P n − n ( t + (cid:15), t ) (cid:15) (A.14)which is exactly the element Q n − n ( t ) of (A.8).In summary we list all the elements of Q ( t )Q nn ( t ) = − [ ν n ( t ) + µ n ( t )] , Q n − n ( t ) = µ n ( t )(A.15)Q n +1 n ( t ) = ν n ( t ) , Q n (cid:48) n ( t ) = 0 , | n (cid:48) − n | > d t φφφ ( t ) = Q ( t ) φφφ ( t ) (A.16)where φφφ ( t ) = ( P ( t ) , · · · , P N ( t )) (cid:124) and Q ( t ) = { Q n (cid:48) n ( t ) } is given by (A.15).It should also be pointed out that the coefficients µ and ν N must be zero, otherwise we require more statesthan [0 , N ] in S . In other words, if µ or ν N are not equalto zero the left side of (A.9) is not equal to one. A2 HOMOEGENEOUS HYPOTHESIS
In this section we show how the hypothesis of spatialand temporal homogeneity are used working in the di-lated sample space.We have already said that if there exists a certain T − generation that can be considered as the stop of theprocess, the upper limit of the dilated space S , denotedby N , can be chosen as N = n T , (A.17)where n T = sup { n : n ∈ S t , ∀ t ∈ [0 , T ] } .Also we considered a dilution of S t into N , i.e. thegeneric dilated space S is equal to N , or N −→ ∞ . Thisassumption is true from S t ⊆ N , no matter how big is n T . In this case the space-time on the process is infiniteand we have an infinitesimal generator on N indepen-dent of the state of the random variable X and time-independent, and the process is space-time homogeneous .However the effect of these hypotheses can be compen-sated with an interaction with the process on N − S t , seethe section A4 and for more details.Essentially we will say that the renormalized distribu-tion defined on N is equivalent to the distribution de-fined on S t . This equivalence is based on the invarianceof evolution equation. In this way both distribution cor-responds to a Markov process, see section A4 . A3 ON THE SPACE-TIME HOMOGENEOUSSOLUTION IN N The evolution equation in the dilated sample space N under the space-time homogeneity is d t φφφ ( t ) = Q φφφ ( t ) , (A.18)with φφφ ( t ) = ( P ( t ) , P ( t ) , · · · ) and identify the matrix Q as Q = − ν µ . . .ν − µ − ν µ . . . ν − µ − ν . . . ν . . . . . . ... ... ... . . . . (A.19)0Under the initial condition φφφ (0) = (1 , , . . . ) (cid:124) , the so-lution of (A.18) with (A.19) is P n ( t ) = e − ( ν + µ ) t (cid:34) ρ n/ I n ( ζ t ) + ρ ( n − / I n +1 ( ζ t )(A.20)+ (1 − ρ ) ρ n ∞ (cid:88) j = n +2 ρ − j/ I j ( ζ t ) (cid:35) . where I n ( x ) is the modified Bessel function [15], ρ = ν/µ and ζ = 2 √ νµ . The first solution of (A.18) appearedin the 1950’s, see [22–25]. A description to obtain thesolution (A.20) is also presented in [14]. A4 MOMENTS OF THE DISTRIBUTION
In this section we calculate the first two cumulants ofthe distribution obtained in (A.20).First of all we demonstrated the existence of all k -moment of the distribution P n ( t ) defined by (cid:104) X k ( t ) (cid:105) = (cid:88) n ∈ N n k P n ( t ) (A.21)with k ∈ N . It is possible to demonstrate that all seriesdefined above converge uniformly ∀ t . To demonstratethis, we define the generatrix function g ( t, z ) = (cid:88) n ∈ N P n ( t ) z n (A.22)for z ∈ R . If it converges, g ( t, z ) is well defined. For ourcase, we know that (cid:88) n ∈ N I n ( x ) −→ [ e x + I ( x )] , (A.23)converges uniformly [15]. This allows us to write thedistribution’s norm and demonstrate that converge uni-formly (cid:88) n ∈ N P n ( t ) −→ . (A.24)Another argument for the general birth-death process,based on the nature of the coefficients { ( ν n , µ n +1 ) : n ∈ N } leads us to the same conclusion [17].For each t , we demonstrated that g ( t, −→ (cid:88) n ∈ N P n ( t ) z n −→ g ( t, z ) (A.25) uniformly for each z ∈ [0 , ∀ ( t, k ) (cid:88) n ∈ N n !( n − k )! P n ( t ) −→ ∂ z k g ( t, z ) | z =1 . (A.26)Finally (cid:104) X k ( t ) (cid:105) can be obtained as a combinationof the series given by (A.26), which converges to ∂ z k g ( t, z ) | z =1 , and this completes the demonstration.We can use two methods for the calculation of the firsttwo moments, given the solution (A.20).The first one is by definition (A.21) and using the iden-tity [15] nI n ( x ) = x I n − ( x ) − I n +1 ( x )] . (A.27)The second one is based on the series (A.21) convergesuniformly. We can derive term by term the series (A.21).Using the evolution equation (A.18) with (A.19), we canobtain a differential equation for the expectation value (cid:104) X ( t ) (cid:105) and (cid:104) X ( t ) (cid:105) d t (cid:104) X ( t ) (cid:105) = ( ν − µ ) + µP ( t ) , (A.28) d t (cid:104) X ( t ) (cid:105) = 2 ν − d t (cid:104) X ( t ) (cid:105) + 2( ν − µ ) (cid:104) X ( t ) (cid:105) . (A.29)For the initial condition P n (0) = δ n then (cid:104) X k (0) (cid:105) = 0and integrate the last two equations (cid:104) X ( t ) (cid:105) = ( ν − µ ) t + µ (cid:90) t P ( τ ) dτ, (A.30) (cid:104) X ( t ) (cid:105) = 2 νt − (cid:104) X ( t ) (cid:105) + 2( ν − µ ) (cid:90) t (cid:104) X ( τ ) (cid:105) dτ. (A.31)The expression (A.30) shows that to determine (cid:104) X ( t ) (cid:105) is sufficient to know the distribution of probability of noblood relationship P ( t ) from (A.20) . Also the expres-sion (A.31) shows that to determine (cid:104) X ( t ) (cid:105) is sufficientto know (cid:104) X ( t ) (cid:105) .This two methods arrives at the same result for (cid:104) X ( t ) (cid:105) .We are interested in computing the first two cumu-lants. The first one is the expectation value and the sec-ond one is defined as function of the first two momentsby (cid:104) [ X ( t ) − (cid:104) X ( t ) (cid:105) ] (cid:105) = (cid:104) X ( t ) (cid:105) − (cid:104) X ( t ) (cid:105) respectively.We have said that this is due to dilation of the samplespace S t (cid:55)−→ N . This implies that distributions definedon N takes smaller values than they should take on S t .1Also we show how (cid:104) X ( t ) (cid:105) is small compared to 2 t +1 ,for any value of ν ≥ µ ≥ ≤ P ( t ) ≤ |(cid:104) X ( t ) (cid:105)| ≤ | ν − µ | t + µ (cid:90) t P ( τ ) dτ ≤ | ν − µ | t + µt. (A.32)This shows that (cid:104) X ( t ) (cid:105) is subordinated to a linear func-tion in t , for all µ, ν , then the expected value of ancestors (cid:104) Y ( t ) (cid:105) = 2 t +1 − (cid:104) X ( t ) (cid:105) it will grow indefinitely with t as2 t +1 . We can consider by ignorance ν = µ , in the senseof not knowing the functional form of the trend of thenumber of ancestors. Although the knowledge of anyparticular trend can be introduced in the gauge trans-formation. And finally without loss of generality we cantake µ = 1 = ν , since the problem of non saturation willbe solved by dilate the distribution to compensate for thedilution of the sample space, as we shall see in the nextsection. This dilation can be seen as a renormalization.In the first place we considered the case where X ( t ) isdefined on N with no renormalization at all. In this way,the renormalization is interpreted as an operation wherethe correct scale of the interaction is retrieved modifying (cid:104) X ( t ) (cid:105) , as if we had solved the problem in the originalsample space S t . A5 ON THE INTERACTION WITH A FICTICIOUSENVIROMENT
We analyzed one of the main concepts: the dilutionof the sample space involves the study of an interactionbetween the initial sample space with the ficticious en-viroment. The introduction of the dilution takes intoaccount the interaction with the lost boundary of S t , un-der the condition of S is large enough to include S t forall t .We can view the process on S t as the result of a processon S which interact with another process defined on thecomplement set R t = S − S t . This point of view can bedescribed in a mathematical precise sense. We define theassociated vector spaces { S, S t , R t } to the sets { S , S t , R t } ,where φφφ is a vector in S , which dim ( S ) = N + 1. If thesample space S t has n t + 1 number of states, we define ϕϕϕ as the first n t + 1 component of φφφ , i.e. ϕϕϕ is vectorof S t . We expressed these vectors in a canonical basis { eee n } n =0 , ··· , N such that φφφ ( t ) = (cid:88) n ∈ S P n ( t ) eee n (A.33) φφφ ( t ) = (cid:88) n ∈ S t P n ( t ) eee n + (cid:88) n ∈ R t P n ( t ) eee n then we have φφφ ( t ) = ϕϕϕ ( t ) + ψψψ ( t ) (A.34)where eee = (1 , , · · · , (cid:124) , eee = (0 , , · · · , (cid:124) and soforth. For construction S t is orthogonal to R t . Thedimensions of these vector spaces are determined by dim ( S ) = N + 1 and dim ( S t ) = n t + 1.From the equation (A.34) we see, roughly speaking,that the process on S t is the result of interaction betweenthe process on S and R t through ϕϕϕ ( t ) = φφφ ( t ) − ψψψ ( t ) (A.35)We can write down the evolution equation (A.16) inthe form d t ϕϕϕ ( t ) = Q ss ( t ) ϕϕϕ ( t ) + Q sr ( t ) ψψψ ( t ) (A.36) d t ψψψ ( t ) = Q rs ( t ) ϕϕϕ ( t ) + Q rr ( t ) ψψψ ( t ) (A.37)where Q ab is the a × b block matrix of Q , for a, b ∈ { s, r } , s = n t + 1 and r = N − n t , explicitly Q = (cid:32) Q ss Q sr Q rs Q rr (cid:33) (A.38)We see that the equations for ϕϕϕ and ψψψ are coupled,hence the interaction character which was noted above.To be more specific, we write a relation of two partialsolutions ψψψ and ϕϕϕψψψ ( t ) = (cid:90) R K rs ( t, t (cid:48) ) ϕϕϕ ( t (cid:48) ) dt (cid:48) (A.39)where (A.39) satisfies (A.37) and K rs ( t, t (cid:48) ) is the kernelof this transformation defined as K rs ( t, t (cid:48) ) = T (cid:26) exp (cid:20) (cid:90) tt (cid:48) Q rr ( τ ) dτ (cid:21)(cid:27) Q rs ( t (cid:48) ) θ ( t − t (cid:48) )(A.40)and θ ( t ) are the unit step distribution and T is the time-ordered operator defined as T [ A ( t ) A ( u )] = (cid:40) A ( t ) A ( u ) : t > u A ( u ) A ( t ) : u > t (A.41)2We can obtain another expression for equation (A.36) d t ϕϕϕ ( t ) = (cid:90) R K ss ( t, t (cid:48) ) ϕϕϕ ( t (cid:48) ) dt (cid:48) (A.42)where K ss is the kernel of the integro-differential equa-tion (A.42) defined as K ss ( t, t (cid:48) ) = Q ss ( t (cid:48) ) δ ( t − t (cid:48) ) + Q sr ( t ) K rs ( t, t (cid:48) ) (A.43)and δ ( t ) is the delta distribution.The same argument allows us to obtain an inverse re-lation of (A.39) ϕϕϕ ( t ) = (cid:90) R K sr ( t, t (cid:48) ) ψψψ ( t (cid:48) ) dt (cid:48) (A.44)simply interchange in (A.40), (A.42) and (A.43) thequantities ϕϕϕ ←→ ψψψ , r ←→ s .The mathematical construction presented here showshow a process in S can be described by the interactionof two sub-processes in S t and R t . Specifically this inter-action can be viewed in the relation (A.39) or its inverse(A.44).In other words, the process on S t can be seen as theinteraction between the processes on S and R t = S − S t ,this interaction emerges from the elements of infinitesi-mal transition probabilities present in K ss ( t, t (cid:48) ) through(A.43) and (A.40). A6 DILATION-DILUTION TRANSFORMATION
The dilution operation of S t into a larger sample space S , can also be understood as an dilation represented bythe substitution rule n t (cid:55)−→ N . In a general sense, every dilution-dilation transformation involves a renormaliza-tion of the distribution obtained above.To illustrate this point let us consider two distributions { p n : n ∈ S } and { P n : n ∈ S } defined over the samplespaces S and S respectively, such that S ⊂ S and bothare normalized in each samples spaces. By definition wehave (cid:88) n ∈ S P n = (cid:88) n ∈ S − S P n + (cid:88) n ∈ S P n = 1 . (A.45)Since all terms of (A.45) are positive, there exist a subsetof S in which the distribution P n , restricted to S , is smallerthan { p n : n ∈ S } . We can see it in this from (cid:88) n ∈ S P n < (cid:88) n ∈ S p n (A.46) We want to describe the process { p n : n ∈ S } throughthe process of { P n : n ∈ S } projecting the distribution P n over S . Since there are values of S for which P n is lessthan p n , that projection should amplify P n to improvethe stated description. This amplification corresponds toa dilation transformation.We had mentioned that the renormalization is per-formed through a linear time-dependent transformation.As we conserve the linearity of the equation (A.16) andsince the original sample space S t is time-dependent, thenorm must be corrected locally. This local transforma-tion is structured as a gauge group , specifically a groupof local dilations .On the other hand, this transformation can be viewedthrough a projection of the distribution P n ( t ) on S t ,which filter a subset of states and renormalizes the distri-bution. We used the concept of conditional probability torelate and motivate the renormalization through a gaugetransformation. To be more clear, we distinguish thegauge transformed distribution p n ( t ) and the projecteddistribution q n ( t ).We used a dilation on a generic sample space S =[0 , N ], for a natural number N .Lets consider the following example. If q n ( t ) is definedas q n ( t ) = P (cid:2) X ( t ) = n | X ( t ) ∈ S t (cid:3) (A.47)this distribution is equal to zero for all state such that n ∈ S − S t . Applying the Bayes identity [17] we have P [ A | B ] P [ B ] = P [ B | A ] P [ A ] (A.48)we can express q n ( t ) as q n ( t ) = Λ( t ) P n ( t ) , (A.49)where Λ( t ) = P (cid:2) X ( t ) ∈ S t | X ( t ) = n ] P (cid:2) X ( t ) ∈ S t (cid:3) . (A.50)The numerator of (A.50) is only equal to 0 or 1 dependingon whether n ∈ S − S t or n ∈ S t , respectively. Thedenominator of (A.50) is P (cid:2) X ( t ) ∈ S t (cid:3) = P (cid:2) X ( t ) ≤ n t (cid:3) and for the mutually exclusive events we have P (cid:2) X ( t ) ∈ S t (cid:3) = (cid:88) n ∈ S t P n ( t ) < q n ( t ) ≥ P n ( t ), for n ∈ S t . The expres-sion (A.49) involves a probability fraction which can beviewed as a projection operator which transforms distri-butions defined on S into distributions defined on S t , byfiltering selected states and redefining the correct nor-malization.The last discussion about the conditional probabilityas a projection operation motivates the following con-struction. Given a process in S described by φφφ ( t ) =( P ( t ) , · · · , P N ( t )) (cid:124) which satisfies a general equation,similar to (A.16) d t φφφ ( t ) = Q ( t ) φφφ ( t ) (A.52)and the dilation transformation g t : P n ( t ) −→ p n ( t ) , (A.53)where p n ( t ) = λ ( t ) P n ( t ) are the components of φφφ (cid:48) ( t ) =( p ( t ) , · · · , p N ( t )) (cid:124) . We showed that the covariance re-quirement is satisfied if g t is a gauge transformation , thisimplies the addition of a covariant derivative. To provethis statement we transform a general evolution equation(A.52) through a local dilation group g t d t φφφ (cid:48) ( t ) − d t λ ( t )[ λ ( t )] − φφφ (cid:48) ( t ) = Q ( t ) φφφ (cid:48) ( t ) . (A.54)In order to keep the shape of (A.16) when we transformedby g t , we need to introduce an affine connection whichcan be implemented through a covariant derivative givenby D t = d t − ω ( t ) (A.55)where ω ( t ) = d t λ ( t )[ λ ( t )] − is the gauge function .In this case, the evolution equation for φφφ (cid:48) ( t ) is D t φφφ (cid:48) ( t ) = Q ( t ) φφφ (cid:48) ( t ) (A.56)which is invariant under the local dilation group.The definition (A.55) allows the distribution P n ( t ) and p n ( t ), related by gauged corresponds to a Markov pro-cess.For the case of this work, the gauge function ω ( t ) isindependent of the t − generation, then λ ( t ) = 2 a t + b , for a , b are real constants. Otherwise, the time dependenceof ω ( t ) introduces a time inhomogeneity in the process. The example used to motivate the definition (A.53)can be extended to describe a rich variety of cases. If wegeneralized q n ( t ) as q n ( t ) = P (cid:2) X ( t ) = n | X ( t ) ∈ S t , Z (cid:3) , (A.57)where Z is an extra condition and eventually space-timedependent. The distribution (A.57) is equal to zero for n ∈ S − S t . Then for the Bayes identity (A.48) appliedto (A.57) q n ( t ) = Λ n ( t ) P n ( t ) (A.58)where Λ n ( t ) = P (cid:2) X ( t ) ∈ S t , Z | X ( t ) = n ] P (cid:2) X ( t ) ∈ S t , Z (cid:3) . (A.59)In a more abstract sense, we consider a partition of S = { B n : n ∈ S } and let us consider the event A n writein the following way: A n = (cid:91) m ∈ S A n ∩ B m . (A.60)Therefore we use the law of total probability P [ A n ] = (cid:88) m ∈ S P [ A n | B m ] P [ B m ] (A.61)for a particular case of the partition B m = { m } and P [ A n ] = q n ( t ), then q n ( t ) = (cid:88) m ∈ S Λ nm ( t ) P m ( t ) (A.62)where Λ nm ( t ) = P [ A n | X ( t ) = m ] . (A.63)But of course, the expression (A.62) contains the exam-ple presented in this work if we take Λ nm ( t ) ≡ λ ( t ) δ nm .The equation (A.62) motivates the natural generalizationof a gauge transformation proposed in (20).4 [1] Derrida, B., Manrubia, S. C., Zanette, D.H., StatisticalProperties of Genealogical Trees. Phys. Rev. Lett. , 9,1987 − PhysicaA −
16 (2000).[3] Derrida, B., Manrubia, S. C., Zanette, D.H., On the ge-nealogy of a population of biparental individual.
J. theor.Biol. , 303 −
315 (2000).[4] Ohno, S., The Malthusian parameter of ascents: Whatprevents the exponential increase of one (cid:112)(cid:112)(cid:112) s ancestors?.
Proc. Natl. Acad. Sci. , 15276 − , 1,311 −
317 (1952).[6] Bennet, J. H., Random mating and sex linkage. J. Theor.Biol. , 1, 28 −
36 (1963).[7] Hardy, G. H., Mendelian Proportions in a Mixed Popu-lation. Sci. N. S. , 49 −
50 (1908).[8] Weinberg, W., ¨U ber den Nachweis der Vererbung beimMenschen . Jahresh. Ver. Vaterl Naturkd. Wrttemb. :369 −
382 (1908), (English translations in Boyer 1963 andJamenson 1977).[9] Begon, M., Townsend, C. R., Harper, J. L.,
Ecology -From Individuals to Ecosystems (Blackwell Publishing,4 th edition 2006).[10] Barucha-Reid, A.T., Elements of the Theory of MarkovProcesses and Their Applications (McGraw-Hill, Seriesin Probability and Statistics 1960).[11] Kijima, M.,
Markov Processes for Stochastic Modeling (Chapman & Hall 1997).[12] Pagel, M., Inferring the historical patterns of biologicalevolution.
Nature , 877 −
844 (1999).[13] Chung, K. L.,
Markov Chains With Stationary TransitionProbabilities (Springer, New York, 2 ◦ edition 1967). [14] Kleinrock, L., Queueing Theory (Volume I, Wiley-Interscience Publication 1975).[15] Abramowitz, M. Stegun, I. A.,
Handbook of mathematicalfunctions (Dover books on mathematics 1972).[16] Feynman, R., Mathematical Formulation of the Quan-tum Theory of Electrodynamics Interaction.
Phys. Rev. , 3, 440 −
457 (1950).[17] Feller, W., An Introduction to Probability Theory andIts Applications. Volume I & Volume 2,
John Wiley &Son Inc. (1968).[18] Serva, M. & Peliti, L., A statistical model of an evolvingpopulation with sexual reproduction,
J. Phys. A: Math.Gen . , 705 −
707 (1991).[19] Liberman, E.; Hauert, C. and Nowak, M. A., Evolution-ary dynamics on graphs.
Nature , 312 −
316 (2005).[20] Heitler, W.,
The Quantum Theory of Radiation (OxfordUniversity Press London, 3 rd edition 1954).[21] Douglas L. T. Rohde, Steve Olson, Joseph T. Chang,Modelling the recent common ancestry of all living hu-mans. Nature , 562 −
566 (2004).[22] N.T.J. Bailey, A continuous time treatment of a simplequeue using generating functions,
J. Roy. Statist. Sot.Ser. B (1954) 288 −
291 (1954).[23] D.G. Champernowne, An elementary method of solutionof the queueing problem with a single server and con-stant parameter,
J. Roy. Statist. Sot. Ser. B − Univ. Michigan , Report M729-IR39, (1953).[25] W. Ledermann & G.E.H. Reuter, Spectral theory for thedifferential equations of simple birth and death processes,