Extremely rare interbreeding events can explain Neanderthal DNA in modern humans
Considering the recent experimental discovery of Green et al that present day non-Africans have 1 to of their nuclear DNA of Neanderthal origin, we propose here a model which is able to quantify the interbreeding events between Africans and Neanderthals at the time they coexisted in the Middle East. The model consists of a solvable system of deterministic ordinary differential equations containing as a stochastic ingredient a realization of the neutral Wright-Fisher drift process. By simulating the stochastic part of the model we are able to apply it to the interbreeding of African and Neanderthal subpopulations and estimate the only parameter of the model, which is the number of individuals per generation exchanged between subpopulations. Our results indicate that the amount of Neanderthal DNA in non-Africans can be explained with maximum probability by the exchange of a single pair of individuals between the subpopulations at each 77 generations, but larger exchange frequencies are also allowed with sizeable probability. The results are compatible with a total interbreeding population of order individuals and with all living humans being descendents of Africans both for mitochondrial DNA and Y chromosome.
The question of whether all of us, living humans, descend exclusively from an anatomically modern African population which completely replaced archaic populations in other continents, or if Africans could have interbred with these local hominids has been the subject of a long lasting and interesting debate. The first of these possibilities, known as Out of Africa model, is based mainly on genetic evidence  further supported by paleontological  and archaeological findings . The latter, known as Multiregional model, on the contrary, has been more supported by morphological studies , but recently it has also been found consistent with genetic data . A third, intermediate possibility, known as assimilation model , suggests that Africans may have interbred with local archaic hominids to a limited extent.
The decision of which model correctly describes the origin of Homo sapiens is obscured by the intricacies of the statistical methods proposed for evaluating the models themselves. Examples of such intricate methods, their conflicting conclusions and subsequent debate are given in [5, 6, 7].
In this paper we will describe by a simple and realistic model the dynamics of two subpopulations – Africans and Neanderthals – interbreeding at a slow rate. In particular, we quantitatively determine the fequency of interbreeding events which are necessary in order that non-African living humans have between 1 and nuclear DNA of Neanderthal origin, according to the discovery of Green et al .
Among other important achievements, the recent seminal paper by Green et al provides the first direct evidence of interbreeding of modern humans with archaic hominids, Neanderthals in this case. By direct evidence we mean having sequenced Neanderthal nuclear DNA and showing that this DNA is more similar to nuclear DNA of living non-Africans than to nuclear DNA of living Africans.
Of course, the findings of Green et al await anxiously for replication by the scientific community. Improvements in the resolution of the genome sequencing, in the comparison with present day individuals and DNA sequencing of other fossils classified as Neanderthals, H. erectus, H. floresiensis and modern humans are mostly welcome.
Based on their findings and on archaeological evidence , it was suggested in  that interbreeding between anatomically modern Africans and Neanderthals might have occurred in the Middle East before expansion of modern Africans into Eurasia, at a time in which both coexisted there. This hypothesis is assumed in this paper, allowing inference of the only parameter in the model, the rate of exchange of individuals between Africans and Neanderthals, and giving some idea on the size of the total population involved in the interbreeding.
The model will be fully explained in the next section, but we anticipate here its main features. Total population size is supposed fixed, but African and Neanderthal subpopulations sizes fluctuate according to the neutral (i.e. Africans and Neanderthals are supposed to have the same fitness) Wright-Fisher model  for two alleles at a single locus. We also assume no biological barriers for interbreeding and no strong hypotheses on the initial composition of the population. Gene flow between subpopulations is implemented by assuming that a fixed number of pairs of individuals per generation is exchanged between them.
The model is characterized by a deterministic component – a system of two linear ordinary differential equations (ODEs) – and a stochastic component – a realization of the Wright-Fisher drift process to be introduced as an external function in the ODEs. The ODEs are exactly solvable, up to definite integrals depending on the stochastic part. The stochastic part can be dealt with by simple simulations.
Assuming a random initial fraction of Africans, our main result is the conditional probability density distribution for the exchange parameter , illustrated in Fig. 1. The condition to be satisifed is that, after interbreeding with Neanderthals, a fraction of 1 to of Neanderthal genes, as suggested by , will be present in the African population. Fig. 1 shows this condition is attained with maximum probability for , i.e. one pair of individuals is exchanged between the two subpopulations every 77 generations. The mean value of is , which corresponds to one pair of individuals exchanged at each 12 generations.
Such conclusions are based on a solvable mathematical model and simple simulations, avoiding statistical in favor of probabilistic methods. Application of probabilistic methods reminiscent of Statistical Mechanics to biological problems has been abundant in the literature of Physics and Mathematics communities, but penetration into Biology and Anthropology has proved more difficult. In particular, both authors of this paper have previously and separately anticipated [11, 12, 13, 14, 15, 16] that evidences based on mitochondrial DNA (mtDNA) could not rule out the possibility of interbreeding among modern humans and other archaic forms. We hope that the direct experimental proof of such interbreeding provided by  can be the occasion for better acceptance of methods such as the ones we will discuss.
While writing the present paper a new report  concerning the interbreeding of modern humans with another archaic hominid group was published. Results have been obtained by studying the fossil nuclear DNA extracted from the finger of a single individual previously known only from its mtDNA . The individual is considered a representative of an archaic group of hominids (Denisovans) different both from moderns and Neanderthals. According to the authors, Denisovan nuclear DNA is present in living Melanesians in a proportion of about . Very few is known about the morphology of Denisovans, as complete fossils belonging to this group are not yet known.
Although we still have no data concerning the size of populations and the duration of coexistence, the model described in this paper might be used to describe the interbreeding between modern humans and Denisovans.
2 The model
Consider a population of constant size equal to individuals. We suppose that the population is divided into two subpopulations we call 1 and 2, generations are non-overlapping and the number of generations is counted from past to future. Reproduction is sexual and diploid. We also suppose that the subpopulations have lived isolated from each other for a long time before they meet. At generation , when subpopulations meet, the total population consists then of two groups, each of which consisting in individuals of a pure race. Starting at this time subpopulations will share a common environment for a long period.
We do not suppose that the numbers and of individuals at generation in each of the two subpopulations are constant, although their sum is. Instead, and are random variables which can be determined by Wright-Fisher rule i.e., any of the individuals of generation independently chooses to belong to subpopulation 1 with probability and to subpopulation 2 with probability . After that, both father and mother of an individual in generation are uniformly randomly chosen among all males and females of generation in the subpopulation he/she has chosen.
With such a reproduction mechanism the numbers and fluctuate as generations pass until one of the subpopulations becomes extinct. This stochastic process is the same as in the simplest version of the neutral, i.e. no selective advantage for any of the alleles, Wright-Fisher model for two alleles at a single locus . The time for extinction is random as well as which of the two subpopulations becomes extinct. If is the initial fraction of individuals of subpopulation 1, then subpopulation 1 will survive with probability and the mean number of generations until extinction is (see ). As the mean number of generations for extinction of one subpopulation scales with , it is reasonable to measure time not in generation units, but in generations divided by . From here on, we will refer to simply as time and we will refer to in a realization of the above stochastic process as the history of the Wright-Fisher drift process.
In the previously described dynamics no mechanism of gene admixture between subpopulations was present and we add it as follows. We assume that at each generation a number of random individuals from subpopulation 1 migrates to subpopulation 2 and vice-versa the same number of random individuals from subpopulation 2 migrates to subpopulation 1. In other words, pairs per generation are exchanged. We strongly underline that is a number of order 1, not of order . Migrants will contribute with their genes for the next generation just like any other individual in their host subpopulation. Their offspring, if any, is considered as normal members of the host subpopulation.
The parameter introduced above may be non-integer and also less than 1. In such cases we interpret it as the average number of pairs of exchanged individuals per generation.
By the hypothesis of isolation between subpopulations for a long time before , we may suppose that in many loci the two subpopulations will have different and characteristic alleles. Therefore, we can assume that there exists a large set of alleles which are exclusive of subpopulation 1 and the same for subpopulation 2. We will refer to these alleles respectively as type 1 and type 2. At any time any individual will be characterized by his/her fractions of type 1 and type 2 alleles. We define then as the mean fraction of type 1 alleles in subpopulation 1 at time and as the mean fraction of type 1 alleles in subpopulation 2 at time . The mean here is due to the fact that individuals in subpopulation 1 in general have different allelic fractions, but is calculated by averaging allelic fractions among all individuals in subpopulation 1, and similarly for . Of course and . Similar quantities might have been defined for type 2 alleles, but they are easily related to and and thus unnecessary.
It is now possible to derive the basic equations relating the mean allelic fractions at generation with the mean allelic fractions at generation . In doing so we will make the assumption that the individuals of subpopulation 1 migrating to subpopulation 2 all have an allelic fraction equal to . The analogous assumption will be made for all the individuals of subpopulation 2 migrating to subpopulation 1.
Of course the above assumption of exchanged individuals all having the mean allelic fractions in their subpopulations is a very strong one and it is not strictly true. Nonetheless, it is indeed a very good approximation if is much smaller than . In fact, is the number of generations between two consecutive exchanges of individuals. As the typical number of generations for genetic homogenization in a population of individuals with diploid reproduction and random mating is , see [19, 20, 21, 22], the condition that is much smaller than makes sure that subpopulations 1 and 2 are both rather homogeneous at the exchange times.
The allelic fraction will be equal to plus the contribution of type 1 alleles from the immigrating individuals of subpopulation 2 and minus the loss of type 1 alleles due to emigration. We remind that these loss and gain terms are both proportional to and inversely proportional to the number of individuals in subpopulation 1. Similar considerations apply to . In symbols:
The above equations, after taking the limit, become a system of linear ODEs
We stress here that we think of as a stochastic function obtained by realizing the Wright-Fisher drift, but Eqs. (1) and (2) still hold if is any description of the history of the size of subpopulation 1, be it stochastic or deterministic. For example, the possibility of individuals in subpopulation 1 being fitter than individuals in subpopulation 2 has been explored, still using (1), in another work .
Eqs. (2) can be exactly solved up to integrals depending on . Although such integrals cannot be calculated in general, the exact solution can be used to give a qualitative view of the behaviour of functions and . It turns out that is a decreasing function, whereas increases. The decrease and increase rates are larger when is large and, despite symmetry in our immigration assumption, gene flow between subpopulations is in general asymmetrical. Such features are shown in appendix A.
Moreover, Eqs. (1) lend themselves to simple and rapid numerical solution for quantitative purposes. In appendix B, we address the question of comparing numerical solutions of Eqs. (1) and direct simulation of all stochastic processes involved. We see that there is good agreement between simulations and numerical solutions of Eqs. (1). In all that follows, unless explicitly stated, we will use results obtained by numerically solving Eqs. (1), because the computer time for numerical solution is much smaller than for simulation.
3 Estimating the exchange parameter
We know that Neanderthals were extinct and, according to , before disappearing they interbred with modern humans. Despite comparisons between nuclear DNA of Neanderthals and living humans having been limited up to now by a sample of only 3 Neanderthals and 5 living humans, the authors of  observed that all three non-Africans in their sample are equally closer to the Neanderthals than the two Africans. They estimate that non-African living humans possess to of their nuclear DNA derived from Neanderthals. Supposing that Africans are subpopulation 1 in our model, this means that the final value of should lie between and in order to comply with their experimental conclusions. We will refer in the following to the interval between 0.96 and 0.99 as the experimental interval for the final value of .
As we do not know the composition of the total population at the time the two subpopulations met, we will take the initial fraction of Africans as a random number. With this hypothesis, the only free parameter is the exchange rate .
As can be seen in Fig. S1 the value of largely influences the final value of . Furthermore, in both Figs. S1 and S2 it can be seen that with or the final values of tend to be too small to be compatible with the experimental interval. We stress that these figures are based only on two realizations of the history and a single value . In order to produce estimates of we must produce a large number of histories with many values of and for any of these simulated histories recursively solve Eqs. (1) in order to determine the associated final value of .
The inset in Fig. 1 is realized by producing 400,000 Wright-Fisher drift histories with random uniformly distributed between 0 and 1. For all these histories we compute the final theoretical value of by solving Eqs. (1) using the three values , and . Therefore, for each of the three values of we have about 200,000 data which allow inference of the probability density for the final value of . The data plotted in the inset of Fig. 1 show that for the probability that the final value of lies in the experimental interval is approximately equal to . For the corresponding probability is approximately of and for it is approximately of . In all three cases the density of the final values of is rather thick, meaning that there is a large probability that the final value of does not lie in the experimental interval.
The above information shows that the experimental data are better explained by values of much smaller than 1. By Fig. 1 we see that the value of which explains with largest probability the experimental data is . In order to produce that plot, we simulated a large number of Wright-Fisher histories with random uniformly extracted between 0 and 0.8 and random values for uniformly distributed between 0 and 2. From these data we selected the histories in which subpopulation 2 was extinct and such that the final theoretical value of lied in the experimental interval. In this way we can empirically determine the probability that the final value of lies in the experimental interval as a function of .
We also see that the probability density for is rather asymmetrical around with values contributing with large probability. This asymmetry is reflected in the fact that the mean value is , much larger than .
A technical detail in producing Fig.1 is that the random values for are chosen with uniform distribution in the interval , avoiding values either close to or inside the experimental interval. Such a choice is related to the assumption of slow rather than rapid interbreeding between Africans and Neanderthals. See aapendix C and Fig. S3 for a more detailed explanation on that choice.
4 Other results
O. Bar-Yosef  compares occupation of the Middle East by Neanderthals and Africans with a long football game. The occupants of the caves of Skuhl and Kafzeh in Israel alternated between Africans and Neanderthals several times over a period of more than 130,000 years. Although the model described before becomes independent of the total population , we may obtain some hints on the size of if we accept the constraint that at least for 130,000 years Neanderthals had not been extinct in the Middle East.
By taking random values for between 0 and 0.8 and between 0 and 2 we obtained a sample of 790 events such that Neanderthals were extinct and lied in the experimental interval. For each of these events we recorded the time it took for extinction of Neanderthals and we found out that the mean extinction time was 0.58. If we take this mean value as the typical value, suppose that one generation is 20 years and equate it to 130,000 years, we get individuals. The whole distribution of extinction times in the above sample is shown in Fig. S4.
In Fig. S3 we plotted the same sample of events in the plane . We see that smaller values of are correlated with smaller values of and also that the events such that lies in the experimental interval are concentrated around the largest values of . The mean value of for the whole sample is 0.64.
Using the same sample we may also explore the values of at the time Neanderthals were extinct, i.e. the fraction of African DNA in the last Neanderthals which interbred with Africans. Fig. S5 shows a histogram of the values for the events in the sample. Observe that typical values of are much larger than the values of , which range from 0.01 to 0.04. This is due to the fact that in most events such that falls within the experimental interval, Africans were the majority of population for most of the time. According to the explanation in appendix A, this implies that, despite symmetry in the number of exchanged individuals, transfer of African alleles to Neanderthals will be larger than transfer of Neanderthal alleles to Africans.
By simulating the complete reproduction and individual exchange process described in appendix B we were able also of empirically determining the conditional probability – the condition being that the fraction of African DNA in Africans is in the experimental interval – that the most recent common ancestors in the population for the maternal (mtDNA) and paternal lineages (Y chromosome) are both African. We ran several simulations with populations of 100 individuals and random values of uniformly distributed between 0.01 and 0.2 and random constrained to be smaller than 0.8. In each simulation we waited until all male individuals had the same paternal ancestor and all female individuals had the same maternal ancestor. We selected those simulations in which subpopulation 1 survived and lied in the experimental interval. Out of 96 simulations satisfying the above criteria, only in 7 of them the survived Y chromosome and mtDNA lineages were not both of ancestors belonging to subpopulation 1. Therefore, according to our interbreeding model, the conditional probability of an African origin of both mtDNA and Y chromosome can be estimated to be of order .
5 Discussion and conclusions
Large samples of mtDNA  and Y chromosomes  in living humans have been sequenced. The small variation among living humans is compatible with a single ancestor woman (mtDNA) and a single ancestor man (Y chromosome) to the whole population, probably both of African origin and living about 100-200 thousand years ago. These facts have been interpreted as proofs of the Out of Africa model, but our interbreeding model is perfectly compatible with them. In fact, conditioned to being in the experimental interval, our model yields a large probability of 93% for African origin of both mtDNA and Y chromosome.
More recently , the whole mtDNA of a few Neanderthal fossils became available. The average number of pairwise differences in mtDNA between a Neanderthal and a living human is significantly larger than the average number of pairwise differences in mtDNA among living humans. This has been considered as a further confirmation of the claim that Neanderthals belong to a separate species, see e.g. , and also for the Out of Africa model.
Before any data on Neanderthal nuclear DNA was available, both authors of this paper had separately anticipated [14, 15, 16, 11, 12, 13] that the above facts are all compatible with anatomically modern Africans and Neanderthals being part of a single interbreeding population at the times they coexisted. Some further details about these claims are given in appendix E.
In the framework of the model proposed in this article we could infer that the 1 to 4% fraction  of Neanderthal DNA in present day non-Africans can be explained with maximum probability by assuming that the African and Neanderthal subpopulations exchanged only 1 pair of individuals in about 77 generations. But the mean value of the exchange parameter in the model corresponds to a larger frequency of about 1 pair of individuals exchanged in about 12 generations.
We also estimated the mean number of generations for Neanderthal extinction in the Middle East to be approximately . Together with the fact that Neanderthals and Africans seem to have coexisted in the Middle East for at least 130,000 years, this allows us to estimate the total population in the model to be of order individuals.
Although Green et al have observed in  gene flow from Neanderthals into Africans, they have not observed the reverse flow. This fact is also compatible both with our results and the fact that living Europeans are as close to Neanderthals as living Asians or Oceanians. The explanation is that the Neanderthal specimens which had their DNA sequenced in  were all excavated in European sites. It seems that only a part of the total Neanderthal population took part in the interbreeding process in Middle East, the other part of the population remaining in Europe. The descendants of these Neanderthals which have never left Europe did not interbreed later with Africans when they came into Europe, or this interbreeding was very small. On the contrary, according to our model, see Fig. S5, we expect to find a larger fraction of African DNA in late Middle East Neanderthal fossils than the 1 to 4% Neanderthal fraction of present non-Africans. Thus, DNA sequencing of one such fossil would be a good test for the present model.
Neanderthals are implicitly considered in this work as a group within the Homo sapiens species and we renounce the strict Out of Africa model for the origin of our species, in which anatomically modern Africans would have replaced without gene flow other hominids in Eurasia. In particular, our model is neutral in the sense that we assign the same fitness to Neanderthals and Africans. Our results show that neither strong sexual isolation between Africans and Neanderthals or else some kind of Neanderthal cognitive or reproductive inferiority, are necessary to explain both their extinction and the small fraction of their DNA in most living humans. In fact, within the assumptions of the model, if two subpopulations coexist in the same territory for a sufficiently long time, only one of them survives. The fact that Neanderthals were the extinct subpopulation is then a random event.
Although we do not intend to back up any kind of superiority for Neanderthals, our neutrality hypothesis is at least supported by recent results [28, 29] by J. Zilhão et al, which claim that Neanderthals in Europe already made use of symbolic thinking before Africans arrived there.
Current knowledge about Denisovans morphology and life style is much less than what we know about Neanderthals. In particular we do not know whether Denisovans lived only in Siberia, where up to now the only known fossils have been found, or elsewhere. Where and when this people made contact with the African ancestors of present day Melanesians is still a mystery. Nevertheless, if such a contact occurred for a sufficiently long time in a small geographical region, then the present model can be straightforwardly applied.
As we now know of our Neanderthal and Denisovan inheritances, it is time to ask whether they were the only hominids that Africans mated. We believe that the future may still uncover lots of surprises when Denisovans will be better studied and nuclear DNA of many more Neanderthal and other hominid fossils will become available.
-  R. L. Cann, M. Stoneking, and A. C. Wilson. Mitochondrial DNA and human evolution. Nature, 325:31, 1987.
-  C. Stringer. Human evolution: Out of Ethiopia. Nature, 423:692–695, 2003.
-  S. Mcbrearty and A. S. Brooks. The revolution that wasn’t: a new interpretation of the origin of modern human behavior. J. Hum. Evol., 39(5):453 – 563, 2000.
-  A. G. Thorne and M. H. Wolpoff. The multiregional evolution of humans, revised paper. Sci. Am., 13(2), 2003.
-  A. R. Templeton. Haplotype trees and modern human origins. Yearb. Phys. Anthropol., 48:33–59, 2005.
-  N. R. Fagundes, N. Ray, M.Beaumont, S. Neuenschwander, F. M. Salzano, and L. Bonatto, S. L. Excoffier. Statistical evaluation of alternative models of human evolution. P. Natl. Acad. Sci. USA, 104(45):17614–17619, 2007.
-  A. R. Templeton. Coherent and incoherent inference in phylogeography and human evolution. P. Natl. Acad. Sci. USA, 107(14):6376–6381;, 2010.
-  R. E. Green, J. Krause, et al. A draft sequence of the Neandertal genome. Science, 328:710–722, 2010.
-  O. Bar-Yosef. Neandertals and modern humans in Western Asia. In T. Akazawa, K. Aoki, and O. Bar-Yosef, editors, The Chronology of the Middle Paleolithic of the Levant, pages 39–56. Plenum, New York, 1999.
-  W.J. Ewens. Mathematical Population Genetics. Biomathematics (Berlin). Springer-Verlag, 1979.
-  A. G. M. Neves and C. H. C. Moreira. Applications of the Galton-Watson process to human DNA evolution and demography. Physica A, 368:132, 2006.
-  A. G. M. Neves and C. H. C. Moreira. The mitochondrial Eve in an exponentially growing population and a critique to the Out of Africa model for human evolution. In R. P. Mondaini and R. Dilão, editors, BIOMAT 2005. World Scientific, 2005.
-  A. G. M. Neves and C. H. C. Moreira. The number of generations between branching events in a Galton-Watson tree and its application to human mitochondrial DNA evolution. In R. P. Mondaini and R. Dilao, editors, BIOMAT 2006. World Scientific, 2006.
-  M. Serva. Lack of self-averaging and family trees. Physica A, 332:387 – 393, 2004.
-  M. Serva. On the genealogy of populations: trees, branches and offspring. J. Stat. Mech.-Theory E., 2005(07):P07011, 2005.
-  M. Serva. Mitochondrial DNA replacement versus nuclear DNA persistence. J. Stat. Mech.-Theory E., 2006(10):P10013, 2006.
-  D. Reich, R. E. Green et al. Genetic history of an archaic hominin group from Denisova cave in Siberia. Nature, 468(7327):1053–1060, 2010.
-  J. Krause, Q. Fu, J. M. Good, B. Viola, M. V. Shunkov, A. P. Derevianko, and S. Pääbo. The complete mitochondrial DNA genome of an unknown hominin from southern Siberia. Nature, 464:894–897, 2010.
-  B. Derrida, S. C. Manrubia, and D. H. Zanette. Statistical properties of genealogical trees. Phys. Rev. Lett., 82(9):1987–1990, Mar 1999.
-  B. Derrida, S. C. Manrubia, and D. H. Zanette. Distribution of repetitions of ancestors in genealogical trees. Physica A, 281(1-4):1 – 16, 2000.
-  B. Derrida, S. C. Manrubia, and D. H. Zanette. On the genealogy of a population of biparental individuals. J. Theor. Biol., 203(3):303 – 315, 2000.
-  J. T. Chang. Recent common ancestors of all present-day individuals. Adv. Appl. Prob., 31(4):1002–1026, 1999.
-  A. G. M. Neves. Interbreeding conditions for explaining Neandertal DNA in living humans: the nonneutral case. Preprint (2011) submitted to BIOMAT 2011, available at http://arxiv.org/PS_cache/arxiv/pdf/1105/1105.6069v1.pdf
-  B. Harder. Did humans and Neandertals battle for control of the Middle East? National Geographic News, March 8 2002.
-  P. A. Underhill, Li Jin, A. A. Lin, S. Q. Mehdi, T. Jenkins, D. Vollrath, R. W. Davis, L. L. Cavalli-Sforza, and P. J. Oefner. Detection of numerous Y chromosome biallelic polymorphisms by denaturing high-performance liquid chromatography. Genome Res., 7:996–1005, 1997.
-  M. Krings, C. Capelli, F. Tschentscher, H. Geisert, S. Meyer, A. von Haeseler, K. Grossschmidt, G. Possnert, M. Paunovic, and S. Pääbo. A view of Neandertal genetic diversity. Nat. Genet., 90(26):144 – 146, 2000.
-  M. Currat and L. Excoffier. Modern humans did not admix with Neanderthals during their range expansion into Europe. PLOS Biol., 2(12):e421, 2004.
-  J. Zilhão, D. E. Angelucci, et al. Symbolic use of marine shells and mineral pigments by Iberian Neandertals. P. Natl. Acad. Sci. USA, 107(3):1023–1028, 2010.
-  K. Wong and J. Zilhão. Did Neandertals Think Like Us? Sci. Am., 302(6):72–75, 2010.
Appendix A Solution and qualitative behaviour of solutions of the model equations
By introducing the auxiliary functions and and taking into account the initial conditions , , we may solve ODEs (2) in the main text of the paper, obtaining
where in Eq. (S2), is given by Eq. (S1). The same path could be followed for the direct solution of the difference equations Eq. (1) in the main text, but formulae corresponding to Eqs. (S1) and (S2) here are more involved and, more importantly, the limit will be appropriate for our further analysis. Of course Eqs. (S1) and (S2) may be trivially used to derive explicit expressions for and , but we think the result is clearer in the form given by Eqs. (S1) and (S2).
In general, is a complicated function obtained by realizing the Wright-Fisher drift. In the limit, it is a solution of the stochastic ODE
where is standard Brownian motion, i.e. and . As a consequence, we cannot explicitly compute the integrals in Eqs. (S1) and (S2). Anyway, Eqs. (S1) and (S2) can be used to give a qualitative description of the solutions to Eq. (2) in the main text and, if necessary, integrals may be easily numerically computed.
As the integrand in the exponent of Eq. (S1) is positive, it shows that the difference between and is positive and steadily decreasing. Moreover, this information, when plugged into Eqs. (2) shows that in fact decreases and increases.
Eq. (S2) on the other hand shows that gene flow from one subpopulation into the other is generally not symmetric. In fact, measures the difference between the fraction of type 1 alleles in subpopulation 2 and type 2 alleles in subpopulation 1. By Eq. (S2), this difference decreases at times in which and increases when . Moreover, it shows that gene flow is more effective at initial times, when values are larger.
Appendix B Checking accuracy of the model and its stochastic simulation
With the purpose of illustrating the qualitative behavior of the solutions of Eqs. (2), see appendix A, we show in Fig. S1 plots of and numerically obtained in the case of two deterministic histories which illustrate typical situations occurring in the Wright-Fisher drift.
It can be seen that all qualitative features of the solutions to Eqs. (2) are present. It should also be noticed in Fig. S1 that the final values of and , i.e. their values at the time of extinction of one of the subpopulations, do depend very much on the history and on the value of .
The final values of and are the most important outputs of the model, because they can be compared with experimental data. As stated above, these values are expected to heavily depend on the particular realization of and on . Therefore, although the qualitative behavior of and , as outlined in appendix A, is quite well-understood, it is necessary to simulate the model by a computer program to obtain quantitative information on their final values.
We first simulate the history . This part begins by choosing a value for the total population size and a value for the initial fraction of individuals in subpopulation 1. Some comments on the choice of are made in appendix C. The choice of is not so relevant if it is large enough so that agreement between the solutions of Eqs. (1) and Eqs. (2) is good. In all results shown we have taken , which produced a good agreement.
Then, individuals in generation independently and randomly choose the subpopulation to which they belong, being the probability of choosing subpopulation 1. The fraction of individuals which expressed the choice for subpopulation 1 produces the value . In general, individuals in generation randomly and independently choose their subpopulation, being the probability of choosing subpopulation 1. This procedure is repeated many times and it generates a realization of the Wright-Fisher drift, i.e. a sequence of values , until attains either the value 0 or the value 1. If a realization of is directly plugged into the difference equations (1), the theoretical values of and can be easily obtained. These theoretical values have been used in producing e.g. results shown in Fig. 1, which represents the core of the paper.
The second part of the program concerns the processes of diploid reproduction and exchange of individuals between subpopulations. In order to obtain the simulated values of and , it is necessary to numerically run the stochastic processes of reproduction and individual exchange. The former consists in the choice of the parents of each individual in the next generation and the latter in the random extraction of individuals to be exchanged between subpopulations. We will explain later the details of them. The purpose of this second part is twofold: it is necessary to check, see Fig. S1, the accuracy of the approximations made in deducing Eqs. (1), and also to obtain information concerning the common ancestors of all individuals in the population in paternal and maternal lineages. Although this information has no relevance for the the final values of and , it is necessary in order to check whether or not the common ancestors of the whole population in paternal and maternal lines belong to the ancestors of the surviving subpopulation.
In the second part of the program, at all time steps we suppose that half the number of individuals in any subpopulation are males and the other half females. The process of individuals exchange is simulated by randomly picking individuals of subpopulation 1 and individuals of subpopulation 2 and exchanging their subpopulation affiliation. In the more interesting case in which is less than 1, we promote the exchange of 1 random individual of each subpopulation each generations.
The reproduction process is simulated as follows: each individual of subpopulation 1 at time makes a random choice of both his/her parents among male and female individuals in subpopulation 1 at time , migrants included. The analogous procedure is followed by individuals in subpopulation 2. At each generation we keep track of the entire genealogy of each of the individuals by counting the number of times each one of the ancestors (individuals of the founding population which lived at time 0, before interbreeding started) appears in his/her genealogical tree . Then we proceed to computing the simulated value of the fractions and . We first consider a single individual at time in subpopulation 1 and we count the number of times the ancestors belonging to subpopulation 1 appear in his genealogy, then we divide this number by , and finally we average this value with respect to all individuals of subpopulation 1 at time . The result is the simulated value of . An analogous calculation produces the simulated value of .
For each male individual at each generation we also keep track of his ancestor by paternal line in generation 0. For the female individuals we do the same for the maternal ancestor in generation 0.
The left graph in Fig. S2 shows the result of one such simulation, in which we compare the theoretical and simulated values for and using the same Wright-Fisher drift history . It should be noted that, although not complete, agreement between simulated and theoretical quantities is good. We remind here that the simulated allelic fractions are subject to statistical fluctuations due to the random processes of exchange of individuals and diploid reproduction.
Indeed, we believe that the randomness in the diploid reproduction process accounts for the largest part of the difference between theoretical and simulated values. In fact, as shown in [19, 20, 21], with diploid reproduction the contribution of each single individual to the gene pool some generations later is highly variable. On the other hand, if is much less than , randomness in the process of exchange of individuals is not so important since at the time of exchanges the individuals in each subpopulation are already highly homogeneous from the point of view of allelic fractions. We have directly checked this fact while producing the data shown in Fig. S2.
It should also be noticed that agreement between theoretical and simulated values is worst for when subpopulation 2 is close to extinction. In this case, in fact, given the small size of subpopulation 2, even a small number of migrants induces large fluctuations in .
The right graph in Fig. S2 shows the average of the simulated values and over 100 simulations with the same history. Notice that the difference between theoretical and average simulated values is accordingly smaller.
Appendix C Why we must exclude values of close to or in the experimental interval when inferring the value of
As remarked in the main text, in producing Fig. 1, we have taken random values of uniformly distributed between 0 and 0.8. The reason why we avoided larger values of is that they are too close to the experimental interval or inside that interval. We now explain why this must be done.
First we observe that if is in the experimental interval, then the final values of and will necessarily also lie in the experimental interval provided that is large enough. The free mating situation, in which subpopulations interact as if there were no differences among their members, is a particular case of this large regime. Free mating, in the infinite population limit, is in fact “described” by Eqs. (2) with an infinite value for . In this case the solution to the equations is straightforward: both and become instantaneously equal to .
The conclusion is that if lies in the experimental interval, then the model would fail to predict any upper bound to , as easy or free mating situations are allowed. Nevertheless, we do not believe that either of these situations were likely to have occurred in reality, since distinct subpopulations coexisted for thousands of years. Therefore, the experimental interval has to be excluded in the choice of .
If we take instead values of outside the experimental interval, but still close to its boundaries, simulations show that both and take very large values, such values tending to infinity as gets closer to the experimental interval. This is illustrated in Fig. S3. With typical values of become comparable to (with the value we used and also with of the order of tens of thousands as it could have been in the real events in Middle East) or larger. As already commented, for such large values of , Eqs. (1) or (2) do not describe accurately the interbreeding process. The reason is that the assumption that all individuals in each subpopulation are genetically homogeneous, necessary to deduce Eqs. (1), fails.
Our choice in Fig. 1 of taking the values of limited to 0.8 is thus a reasonable consequence of the mathematical characteristics of the model. It is also a reasonable choice from a historical point of view, because we are assuming that the Neanderthal subpopulation was comparable to the African one; it might be smaller, but not extremely smaller, compared with the African one.
Appendix D Other results
The main result of our paper is the probability density distribution for the interbreeding rate shown in Fig. 1. Other interesting resuts are illustrated here. In order to produce them we simulated a large number of Wright-Fisher drift histories with random constrained to be smaller than 0.8 and random smaller than 2. For each history we numerically obtained the values of ad by iterating Eq. (1). We obtained then a set of 790 events such that the final value of lies in the experimental interval.
Fig. S3 was produced with the above mentioned sample. With the same sample we may also study the questions of extinction times, see Fig. S4, and the final values of , i.e. how much African DNA was transmitted to Neanderthals before they were extinct in the Middle East, see Fig. S5. Comments on these figures were made in the main text.
Appendix E Mitochondrial DNA and Y chromosome
Mitochondrial DNA and Y chromosome are both inherited in a haploid way. Furthermore mtDNA is not subject to recombination and recombination seems to be negligible for the Y chromosome. It is also believed that large portions of both are selectively neutral. These facts allow an easier mathematical treatment of their statistical properties. From the experimental point of view, large samples of mtDNA  and Y chromosomes  in living humans have been sequenced. The small variation among living humans is compatible with a single ancestor woman (mtDNA) and a single ancestor man (Y chromosome), probably both of African origin and living about 100-200 thousand years ago. These facts have been interpreted as proofs of the Out of Africa model.
More recently , the whole mtDNA of a few Neanderthal fossils became available. The average number of pairwise differences in mtDNA between a Neanderthal and a living human is significantly larger than the average number of pairwise differences in mtDNA among living humans. This has been considered as a further confirmation of the claim that Neanderthals belong to a separate species, see e.g. , and also for the Out of Africa model.
Both authors of this paper have separately claimed that the above facts are all compatible with anatomically modern Africans and Neanderthals being part of a single interbreeding population at the times they coexisted. In , using Kingman’s coalescence, it was shown that the probability distribution of genealogical distances in a population of fixed size and haploid reproduction is random even in the limit when the population size is infinite. The random distribution typically allows large genealogical distances among subpopulations. In  another important fact was statistically described: in a population of fixed size and haploid reproduction one of the two main subpopulations will become extinct at random times with exponential distribution. When such an extinction occurs, average genealogical distances among individuals in the population have a sudden drop. Finally, in , it was shown that mtDNA may be completely replaced in a population by the mtDNA of another neighbor population, whereas some finite fraction of nuclear DNA persists.
These facts imply that the large genealogical distances between living humans and Neanderthals, as seen in mtDNA, are not uncommon in an interbreeding population. On the contrary, they turn out to be very likely if the correct statistics is used. Furthermore, these facts imply that these distances may have been much larger at the time of Neanderthals extinction than they are nowadays. They also imply that extinction of Neanderthals’ mtDNA is compatible with the survival of their nuclear DNA.
Exactly the same reasoning can be applied to the mitochondrial and nuclear DNAs of the fossil bones found in Siberia [18, 17], later described as the new population of Denisovans. The fact that Denisovans differ significantly both from Neanderthals and living humans in their mtDNA  does not imply that they could not interbreed with either of them. Indeed, nuclear DNA proved  that they have interbred at least with some anatomically modern populations.
In [11, 12] the authors examined the question of survival of mtDNA and Y chromosome lineages in a population subject to exponential stochastic growth (supercritical Galton-Watson branching process). It was shown that exponential growth is compatible with the survival of a single mtDNA or Y chromosome lineage only if the growth rate is in a narrow interval. Thus, even if Neanderthals and anatomically modern Africans belonged to the same interbreeding population and even if this population was allowed to grow exponentially with a small rate, the more probable outcome would still be all humans being descendants either of a single woman (mtDNA) or a single man (Y chromosome).
In , the number of generations between successive branchings in the Galton-Watson process was computed. It was found that in the slightly supercritical regime, in which the survival of a single lineage is expected, trees typically have very long branches of the size of the whole tree along with shorter branches of all sizes. Thus, trees are qualitatively similar to those of the coalescent model and, as a consequence, the phenomenon of sudden drops in genealogical distances, described in  is also present in this model.