## Abstract

Two sequentially Markov coalescent models (SMC and SMC′) are available as tractable approximations to the ancestral recombination graph (ARG). We present a Markov process describing coalescence at two fixed points along a pair of sequences evolving under the SMC′. Using our Markov process, we derive a number of new quantities related to the pairwise SMC′, thereby analytically quantifying for the first time the similarity between the SMC′ and the ARG. We use our process to show that the joint distribution of pairwise coalescence times at recombination sites under the SMC′ is the same as it is marginally under the ARG, which demonstrates that the SMC′ is, in a particular well-defined, intuitive sense, the most appropriate first-order sequentially Markov approximation to the ARG. Finally, we use these results to show that population size estimates under the pairwise SMC are asymptotically biased, while under the pairwise SMC′ they are approximately asymptotically unbiased.

- sequentially Markov coalescent
- ancestral recombination graph
- consistency
- ergodicity
- Markov approximation

OF the many models of genetic variation in the field of population genetics, few have as much relevance in the era of genomics as the ancestral recombination graph (ARG). The ancestral recombination graph models patterns of ancestry and genetic variation within sequences experiencing recombination under neutral conditions (Hudson 1991; Griffiths and Marjoram 1997). Under the formulation of Griffiths and Marjoram (1997), lineages recombine apart and coalesce together back in time to produce a graph structure describing the ancestral genealogy at each point along a continuous chromosome. While only a few simple rules govern the process, many aspects of the model are analytically intractable.

Wiuf and Hein (1999) provided a formulation of the ARG that proceeds across the chromosome (rather than back in time), producing the genealogy at each point sequentially. As with the back-in-time formulation of the ARG, at each point along the chromosome there is a local genealogy describing the ancestry of the sample at that point, and changes in the genealogy occur at recombination sites. In this sequential formulation of the ARG, a new lineage is produced wherever an ancestral recombination event is encountered along the chromosome. To produce a new genealogy at the recombination site, the new lineage is coalesced to the ARG representing the ancestry of all previous points along the chromosome. This dependence on all previous points makes the process non-Markovian along the chromosome and (together with a large state space) makes calculations often intractable.

Approximations to the ARG have been suggested with the goal of modeling coalescence with recombination in a way that is analytically tractable. McVean and Cardin (2005) introduced the *sequentially Markov coalescent* (SMC). The original formulation of the SMC was sequential, generating genealogies along the chromosome such that each new genealogy depends only on the previous genealogy. Like the ARG, the SMC has both a back-in-time formulation and a sequential formulation. The back-in-time formulation of the SMC is equivalent to that of the ARG except that coalescence is allowed only between lineages containing overlapping ancestral material. As a consequence, in the sequential formulation of the pairwise ( chromosomes) SMC, each recombination event produces a new pairwise coalescence time.

Marjoram and Wall (2006) introduced a slight modification to the SMC, termed the SMC′, which retains the Markov behavior along the chromosome but models additional coalescence events that make it a closer approximation to the ARG. Specifically, in the back-in-time formulation of the SMC′, coalescence is allowed between lineages containing either overlapping *or adjacent* ancestral material. In the sequential formulation of the pairwise SMC′, this means that not every recombination event necessarily produces a new pairwise coalescence time, since two lineages created by a recombination event can coalesce back together. Figure 1 shows the transitions that are permitted under the back-in-time and sequential formulations of the pairwise ARG, SMC, and SMC′. The sequentially Markov coalescent models have been used in many recently introduced population-genetic, model-based inference procedures, including the pairwise SMC (PSMC) model (Li and Durbin 2011), the multiple SMC (MSMC) model (Schiffels and Durbin 2014), diCal (Sheehan *et al.* 2013), coalHMM (Hobolth *et al.* 2007; Dutheil *et al.* 2009), and ARGWeaver (Rasmussen *et al.* 2014), and in a study of past human demography based on tracts of identity by state (Harris and Nielsen 2013).

The SMC′ was shown by simulation to produce measurements of linkage disequilibrium more similar to the ARG than those produced by the SMC (Marjoram and Wall 2006) and was hence used as the preferred model by some recent studies (Harris and Nielsen 2013; Schiffels and Durbin 2014; Zheng *et al.* 2014). Additionally, a number of recent studies have explored the theoretical properties of the SMC′ (Eriksson *et al.* 2009; Harris and Nielsen 2013; Carmi *et al.* 2014; Schiffels and Durbin 2014; Zheng *et al.* 2014). However, few direct comparisons between the SMC′ and the ARG have been made, and a number of open questions remain. Here, we show how the joint distribution of pairwise coalescence times at two fixed points along a chromosome evolving under the SMC′ can be described by a continuous-time Markov chain. Through analysis of this Markov chain, we calculate many statistical properties of the pairwise SMC′ and compare them to those of the ARG and the SMC. Specifically, for each model of coalescence with recombination, we compare the following: the joint density (*Joint probability density functions*), the conditional density (*Conditional distribution of coalescence times*), and the covariance between and , which we show to be equal to the probability that and are the same (*Covariance of coalescence times*). These quantities are readily related to measures of linkage disequilibrium in real sequence data.

Using our two-locus Markov process for the two-locus, pairwise SMC′, we also show that the joint distribution of coalescence times immediately to the left and right of a recombination event is the same under the SMC′ and ARG. This allows us to calculate the asymptotic bias of the pairwise SMC- and SMC′-based population-size estimators, which we confirm by simulation. We show that the SMC′ estimator is approximately asymptotically unbiased.

## Results

### Two-locus Markov chain model for the SMC and SMC′

Here, we present back-in-time Markov processes for the two-locus SMC and SMC′. Previous work has developed analogous two-locus, back-in-time Markov processes for the ARG. Kaplan and Hudson (1985) first described how the process of generating coalescence times at two linked loci modeled by the ARG can be represented as a continuous-time Markov chain, with coalescence and recombination events causing transitions between states. Simonsen and Churchill (1997) explored this process further for the case where the sample size is and derived for the ARG many of the quantities we compare against the SMC′ in this article. Subsequent work has extended this approach to study two-locus coalescence distributions in the presence of population structure (Eriksson and Mehlig 2004) and recurrent bottlenecks (Schaper *et al.* 2012) and to study species-tree concordance at linked loci (Slatkin and Pollack 2006) and coalescence histories at one locus conditional on the history at a nearby locus (Hobolth and Jensen 2014).

We begin by presenting the simpler SMC model, which provides context for the more complex SMC′ model. If time is scaled such that the rate of coalescence is 1 and the total rate of recombination between the two linked loci is , then the two-locus ancestral process under the SMC is the process depicted in Figure 2. The process starts in state with two lineages, each containing linked copies of the two loci. From , the process transitions with rate *ρ* to state , in which one of the two chromosomes has experienced a recombination event, or with rate 1 to state , an absorbing state in which both loci have coalesced. Under the SMC, lineages can coalesce only if they contain overlapping ancestral material, so after entering , the process cannot return to the fully linked state , and each locus coalesces independently with rate 1 from that time onward. Thus, under the SMC, the joint distribution of coalescence times at two loci is that of (1)where is the amount of time to leave , indicates whether the first event is a recombination event, and are the exponential waiting times until coalescence after the first recombination event. All of these random variables are independent in the SMC model, so it is straightforward to calculate many of the quantities we compare in this article, using this representation.

The defining rule of the SMC′ model of coalescence with recombination is that only ancestral lineages containing overlapping or contiguous ancestral material can coalesce (Marjoram and Wall 2006). The back-in-time process of coalescence at two fixed loci under this model is the continuous-time Markov chain shown in Figure 3. Under the SMC′, it is necessary to model the number of recombination events that have occurred between the two loci at each point in time. To see that this is the case, consider the state in Figure 3. In this state, two recombination events have occurred between the focal loci, and neither focal locus has coalesced. Because lineages can coalesce only to lineages containing overlapping or adjacent ancestral material, two particular coalescence events would need to occur before the process returns to state , regardless of the placement of the recombination events on the two chromosomes.

The SMC′ two-locus Markov process also features an additional state **I**, which is entered when some portion of the chromosome between the focal loci coalesces before either of the focal loci. Upon entering **I** it becomes impossible for the process to reenter the initial, fully linked state (), so the remaining times until coalescence at the focal loci become independent exponential random variables with mean 1. If is the state in which neither focal locus has coalesced and *i* recombination events have occurred between the focal loci, the transition rate into **I** is . This is due to the fact that each recombination event after the first produces an additional pair of lineages that can coalesce to take the process to **I**. For each state , , the number of lineages that can coalesce to take the process to is *i*, and the rate of transitioning to through recombination is *ρ*. Transitions to and occur at rate 1 whenever the process is in state , . Following Eriksson and Mehlig (2004), we disregard any information about linkage between the two loci after one locus has coalesced, since the rate of coalescence at the uncoalesced locus is 1 regardless of the state of linkage with the coalesced locus.

For comparison, an analogous two-locus continuous-time Markov chain for the ARG is presented in Supporting Information, Figure S1. An equivalent process was studied by Simonsen and Churchill (1997) and others. In this model, state is reached when the first event is a recombination event, and state is reached only after a subsequent recombination event occurs on the ancestral lineage that did not experience the first recombination event, making all ancestral copies of the two loci unlinked.

### Joint probability density functions

Considering the SMC′ model above, let represent the probability that the two-locus ancestral coalescent process is in state at time *t*, and let represent the probability that the process is in any state , , or state **I**, at time *t*. The joint density of coalescence times at the two focal loci is then (2)since is the rate of entering state at time *t*, and is the rate of entering either or at time *t*. The joint density for the ARG and SMC is analogously defined, with representing and under the ARG and under the SMC. For the ARG and the SMC, the number of states is finite and and can be solved using matrix exponentiation. For the SMC′, there is an infinite number of states, representing the possibility of an infinite number of recombination events occurring between the two focal loci. In the *Appendix*, we derive closed-form expressions for , , and . The main idea in these derivations is to recognize the structure of the SMC′ in Figure 3 as a birth–death process with killing. In this formulation the states are , a birth corresponds to a recombination event (and the birth rate is constant), a death corresponds to a coalescence event (and the death rate is linear), and killing corresponds to leaving the states.

Figure 4 compares the joint coalescence time distributions under the SMC and SMC′, displaying the differences of these joint distributions with the joint distribution of the ARG. The SMC′ provides a much better fit to the ARG for the range of recombination rates compared. Both the SMC and the SMC′ underestimate the density of outcomes where , but this underestimation is substantially less under the SMC′.

To summarize the difference between the joint distributions more succinctly, we calculated the total variation distance between the SMC and the ARG and between the SMC′ and the ARG across a range of recombination rates. The total variation distance between the SMC and the ARG is defined as (3)where and are the joint densities defined under the SMC and ARG, respectively. The total variation distance between the SMC′ and the ARG is similarly defined. Figure 5 shows the total variation distance from the ARG for the SMC and SMC′ over a range of recombination rates. Total variation distances were calculated numerically. For both the SMC and SMC′, the total variation distance was maximized at some intermediate recombination rate, for the SMC and for the SMC′.

### Conditional distribution of coalescence times

In this section we consider the distribution of coalescence times at one locus given the coalescence time at the other. The conditional density of given , , can be calculated by dividing the joint density by the marginal distribution of coalescence times at the left locus:

(4)Hobolth and Jensen (2014) introduced a framework for modeling the distribution of given , using a time-inhomogeneous continuous-time Markov chain. [Note that the model called SMC′ in Hobolth and Jensen (2014) is an SMC′-like model of two loci that is not based on the continuous-chromosome SMC′. It is different from the SMC′ model we consider here.] This framework can be extended to the SMC′, producing the continuous-time Markov chain shown in Figure S2. Figure 6 compares the conditional density of coalescence times at the right locus conditioned upon the coalescence times at the left locus for different values of and *ρ*.

We note that recently it was proposed that the mutation rate could be estimated by simulation-based calibration of the increase in mean heterozygosity when moving away from a site of known, low heterozygosity (Lipson *et al.* 2015). Our expressions for the conditional distribution of coalescence times could provide theoretical expectations for such a statistic.

### Covariance of coalescence times

In the two-locus, back-in-time Markov processes for the SMC, SMC′, and ARG, and are equal when the state is entered through rather than or . For the ARG, Simonsen and Churchill (1997) showed that the probability that is equal to is (5)Under the SMC (McVean and Cardin 2005), (6)Eriksson *et al.* (2009) used the sequential formulation of the SMC′ to show that (7)where is the exponential rate of encountering a change in coalescence time when the local coalescence time is *t* and is the incomplete gamma function.

For the ARG and SMC, the covariance is equal to . Eriksson *et al.* (2009) showed by simulation that this is also true of the SMC′. Here we present a short proof that this is the case for any two-locus model of coalescence where the marginal distribution of coalescence times is exponential with rate 1.

The expectation can be derived using the fact that (8)The final equality in (8) follows from the fact that has an exponential distribution with rate 1 when . Therefore and

(9)This result holds in other situations with exponential coalescence times, for example in the context of the population-divergence model considered by Eriksson *et al.* (2009) (in which case the marginal distribution is exponential plus a constant) and for the various covariances used by McVean (2002) to calculate , the approximation to the linkage disequilibrium measure .

It is interesting to consider when *ρ* is small. For the ARG, consideration of (5) shows that . Likewise, for the SMC, (6) shows that .

For the SMC′, the integral representation of in (7) allows for the calculation of this quantity as a first-order expansion in *ρ*: (10)Thus, [or ] is the same up to order under the ARG and SMC′.

### Coalescence times at recombination sites

In this section, we show that the joint distribution of coalescence times on either side of a recombination event is the same under the SMC′ and marginally under the ARG, and we derive this distribution. Consider the continuous-time Markov chains representing the two-locus SMC′ and ARG models (Figure 3 and Figure S1, respectively) in the limit of and conditioning on the first event being a recombination event. These processes represent the joint distribution of coalescence times on either side of a recombination event under the ARG and SMC′. In both of these processes, the waiting time until the first event, conditional on that event being a recombination event, has an exponential distribution with rate , which converges to 1 as After that first recombination event, the rate of all additional recombination events converges to zero in the limit, so all of the remaining events must be coalescence events, each of which occurs with rate 1. Under the ARG and the SMC′, the coalescence events that are possible from state are the same. Thus, the joint distribution of coalescence times at recombination sites is the same under the SMC′ and the ARG.

Figure 7A shows the two-locus continuous-time Markov chain representing this conditional process. This Markov chain starts in a special initial state , out of which the first event is always a recombination event, which happens with rate 1, as described above. In previous sections, we used and to represent the coalescence times at two loci some fixed distance apart. To avoid confusion, in this section we use *S* and *T* to represent the coalescence times on the left and right sides of a recombination event, respectively.

Recombination events are visible in sequence data only if they change the local coalescence time. Thus, it is of special interest to condition on in the above model to derive the joint distribution of coalescence times on either side of a change in coalescence times under the ARG and SMC′. Conditioning on , the transition out of must be into either or . These transitions occur with conditional rate , since the total rate of leaving is 3 in the unconditional model, and two of the ways of leaving result in the coalescence times being different.

The continuous-time Markov chain representing coalescence times on either side of a change in coalescence times (*i.e.*, at recombination sites where ) is shown in Figure 7B. Under this model, the joint distribution of *S* and *T* is that of (11)where , , , , and the random variables are independently distributed.

Under the SMC, the continuous-time Markov chain representing the joint distribution of coalescence times at recombination sites is equivalent to the model in Figure 7B with the transition rates from to and equal to 1 instead of . Under this model for the SMC, the joint distribution of coalescence times on either side of a recombination event is that of (12)where , , and are mutually independent exponential random variables with rate 1.

In File S1, we use these Markov processes to derive the joint, marginal, and conditional distributions of coalescence times at recombination sites under the ARG, SMC′, and SMC. These calculations confirm previous derivations of Carmi *et al.* (2014) for the SMC′ and Li and Durbin (2011) for the SMC.

### SMC′ as the canonical first-order Markov approximation to ARG

Under the sequential formulation of the continuous-chromosome ARG, SMC, and SMC′ models, the infinitesimal probability of a recombination event occurring in the interval given the coalescence time *s* at *x* is . This fact, together with the fact that the joint distribution of coalescence times at recombination sites is the same under the ARG and SMC′ (whether or not the coalescence time changes), implies that the conditional distribution of coalescence times at point given the coalescence time at point *x* is the same under the SMC′ and ARG.

This demonstrates that the pairwise SMC′ is the canonical first-order Markov approximation for the pairwise ARG. Given an infinite-order Markov chain , where the distribution of each depends on all previous , , the canonical *k*-order Markov approximation to is the Markov chain satisfying (13)That is, the transition probabilities under the *k*-order canonical Markov approximation are equal to the transition probabilities conditional on the previous *k* states under the infinite-order chain. See Schwarz (1976), Fernández and Galves (2002), and Gallo *et al.* (2013) for examples of mathematical studies of canonical Markov approximations of infinite-order Markov chains.

Here we informally extend the terminology of canonical Markov approximations to continuous processes. The SMC′ is the canonical first-order Markov approximation to the ARG because the distribution of coalescence times at conditional on the coalescence time at *x* is the same under the ARG (an infinite-order, sequentially non-Markovian continuous process) and the SMC′ (a first-order sequentially Markov continuous process). In this sense, the SMC′ is the most natural first-order sequentially Markov approximation to the ARG.

## Asymptotic bias of the population-size estimators under SMC and SMC′

Given the joint density of pairwise coalescence times at recombination sites under the ARG, it is possible to determine the asymptotic bias of maximum-likelihood population size estimators derived from the pairwise SMC and SMC′ likelihood functions. These likelihood functions give the probability of observing a sequence of pairwise coalescence times and corresponding segment lengths across a chromosome under the SMC and SMC′ models. Related likelihood functions (allowing for variable historical population size) are implicitly maximized in the PSMC and MSMC inference procedures (Li and Durbin 2011; Schiffels and Durbin 2014, respectively). These inference procedures are hidden Markov model (HMM) methods in which the local coalescence times (or genealogies) and segment lengths are hidden states inferred from sequence data.

Here, we consider the estimators that would be obtained if the hidden states in these models were actually observable (see also Kim *et al.* 2015). We are motivated by the fact that any biases of the estimators we investigate are likely to be inherent in the full HMM-based inference procedures, since these biases would be present even with perfect knowledge of an infinite number of coalescence times. Furthermore, by analyzing estimators derived from the hidden coalescence states, we isolate the bias that is due to choice of coalescent algorithm (SMC *vs.* SMC′) from the bias due to the mutation model or discretization of the continuous hidden states in a full HMM approach to inference.

To investigate the asymptotic properties of these estimators, we assume that data are generated under the ARG, such that at a fixed point the distribution of pairwise coalescence times is exponential with rate = 1 and an ancestral segment of length *l* recombines back in time at rate . Segment lengths are measured in units of the true scaled recombination parameter *ρ*. Data generated under this model can be represented as a sequence of pairwise coalescence times and corresponding segment lengths: .

We are interested in estimating a single relative population size *η* (defined relative to the true population size, *N*). If the data are modeled by the SMC or SMC′, the likelihood of a particular value of *η* is (14)where is the transition function and is the rate of encountering the end of a segment given *t*, with both quantities pertaining to the sequentially Markov coalescent model being used to calculate the likelihood.

In the *Appendix*, we show that if the SMC is used, the maximum-likelihood estimate of *η* converges to ∼0.95 as the chromosome gets infinitely long. If the SMC′ is used, the estimate is approximately unbiased in the same limit. If the data are reduced to just the segment ages, the likelihood equation is (15)Using this reduced likelihood, the asymptotic maximum-likelihood estimate is asymptotically unbiased under the SMC′. Under the SMC, the reduced likelihood and the full likelihood produce the same maximum-likelihood estimate (see *Appendix*).

We confirm the asymptotic bias of the SMC estimator and the apparent lack of asymptotic bias of the SMC′ estimators by simulation. Figure 8 shows 100 simulated estimates calculated using the SMC, SMC′, and reduced SMC′ likelihood functions. Each estimate was calculated using 100 independent pairs of chromosomes simulated under the ARG, with each chromosome of total length , where *N* is the diploid size and *r* is the per-generation probability of recombination. Likelihood functions were multiplied across independent pairs of chromosomes, and the same set of simulations was used to produce the estimates for all three likelihood functions.

## Discussion

We have presented a continuous-time Markov chain that describes the pairwise coalescence times at two fixed loci evolving under the SMC′ model of coalescence with recombination. We analyzed this Markov chain to derive the joint distribution of coalescence times at the two loci and the conditional distribution of coalescence times at one locus given the coalescence time at the other. We compared these distributions to those of the ARG and SMC models and found that the difference between the ARG and the SMC′ was much less than the difference between the ARG and the SMC.

We showed that the conditional distribution of coalescence times at point given the coalescence time at *x* is the same under the ARG and SMC′. This implies that the SMC′ is the canonical first-order approximation to the pairwise ARG. However, this correspondence is true only of the continuous-chromosome models. If instead the ARG is a model of the genealogies at a sequence of discrete loci, then the first-order canonical Markov approximation is the Markov approximation obtained by modeling a conditional ARG between every successive pair of loci. This model was studied by Hobolth and Jensen (2014), who referred to the model as a “natural” Markov approximation to the ARG. Conceptually similar sequentially Markov coalescent models have been used in the so-called “coalescent hidden Markov model” family of methods (Hobolth *et al.* 2007; Dutheil *et al.* 2009; Mailund *et al.* 2011).

Chen *et al.* (2009) presented a method of simulating data under higher-order sequentially Markov approximations to the ARG, where the ARG of some number of preceding loci is retained in the process of generating the marginal genealogy at a given locus. They showed by simulation that higher-order approximations generate times until most recent common ancestry that are more consistent with the ARG than do lower-order approximations, but little theoretical work on these higher-order Markov approximations has been done.

The two-locus Markov chains analyzed in this article assume a single well-mixed population, but natural populations often have more complex demographic histories, featuring, for example, variable historical population sizes, migration between subpopulations, and/or past divergence from other populations. The theoretical properties of the sequential, across-the-chromosome formulations of the pairwise SMC and SMC′ with variable population sizes have been studied previously (Li and Durbin 2011; Schiffels and Durbin 2014). Eriksson *et al.* (2009) used simulation to study two-locus properties of the SMC′ with population bottlenecks, migration between subpopulations, and divergence between populations. They found that the SMC′ generally performs well in these contexts. The two-locus Markov chains we study here could be extended to include these features (as was done for the ARG by Lessard and Wakeley 2003 and Eriksson and Mehlig 2004), which would provide a framework for calculating exact quantities for the two-locus SMC and SMC′ in the context of structured populations. We leave this for future work.

We calculated the asymptotic bias of a population size estimator under the pairwise SMC to be ∼95% of the true population size. This is not a very large bias, but given the continued use of the SMC model in population-genomic inference methods (Palamara *et al.* 2012; Sheehan *et al.* 2013; Rasmussen *et al.* 2014), there is an apparent need to reexamine the consequences of using the simpler SMC model instead of the slightly more complicated SMC′ model. For example, it will be important to consider whether including the possibility of varying population sizes will increase or decrease asymptotic bias. In this context, using the SMC as a basis for a likelihood function may also bias the estimates of the magnitude and timing of population size changes, since the longer segments produced by the ARG will seem younger when they are modeled under the SMC.

Depending on the particular application, it may sometimes be mathematically difficult to employ the SMC′ instead of the SMC. Nevertheless, the SMC′ is the model underlying two recently introduced population-genetic inference methods: the MSMC method of Schiffels and Durbin (2014) (which simplifies to a PSMC′ inference procedure when the number of haplotypes is two) and a procedure based on the distribution of distances between heterozygous bases, introduced by Harris and Nielsen (2013). In each case it was acknowledged that the SMC′ provided more accurate results than the SMC.

From the arguments that led to the development of the continuous-time Markov chains representing the joint distribution of coalescence times at recombination sites (Figure 7), it seems that the joint distribution of coalescence times on either side of a recombination event will be the same under a variety of demographic scenarios. If one were to allow the historical population size to vary, the waiting time until the conditioned-upon recombination event would still be the same under the SMC′ and ARG, and the remaining coalescence events would also be distributed identically. Likewise, when there is population substructure with migration between subpopulations, the distribution of events occurring at recombination sites should be the same under the SMC′ and ARG. Finally, when there are more than two haplotypes sampled, it seems that the joint distributions of genealogies on either side of a recombination event would be the same between the SMC′ and the ARG marginally. These ideas need to be properly explored in future studies, but they suggest that asymptotic bias due to using the SMC′ in place of the ARG will be minimal under a variety of demographic scenarios.

## Acknowledgments

We thank Erik van Doorn and Søren Asmussen for identifying the correspondence to birth–death models with killing (see *Appendix*). We are grateful to John Wakeley and Paul Marjoram and two anonymous reviewers for comments that helped improve this article. S.C. thanks the Human Frontier Science Program for financial support.

## Appendix

### Derivation of Joint Density of Pairwise Coalescence Times at Two Loci

To calculate the joint density of coalescence times, it is necessary to calculate , the probability that the SMC′ two-locus Markov process (Figure 3) is in state at time *t*, and , the probability that the SMC′ process is in state **I** at time *t*. To solve for , one can use the forward Kolmogorov equation (for ) (A1)Through substitution, the solution to (A1) can be shown to be (A2)To find , we note that it is equal to (see Equation 2). In turn, (A3)where is the marginal distribution of coalescence times at the first (or second) locus and is the probability of no change in coalescence times given the coalescence time *t* at the first locus. Here is the exponential rate of encountering a change in coalescence time along the chromosome given that the local coalescence time is *t* (Eriksson *et al.* 2009; Carmi *et al.* 2014). Thus is given by (A4)This completes the solution of . Using Figure 3, (A5)where is the probability that the process is in state **I** at time *t*. Using (A2) and (A4), we get (A6)Next, satisfies the forward Kolmogorov equation (A7)the solution to which is (A8)Here, is the incomplete gamma function.

Together (A4), (A5), (A6), and (A8) give the joint distribution (2) for the SMC′. For the ARG and SMC, the quantities analogous to and for these models can be obtained by exponentiating the rate matrices implicit in Figure 2 and Figure S1. For the SMC, the joint distribution can also be derived using the representation (1).

The walk on the states constitutes a birth–death process with killing, where birth events correspond to additional recombination events taking the process from to ; death events correspond to coalescence events that take the process from to ; and killing events, which take the process to an absorbing state, here correspond to coalescence events that take the process to , , or **I**. Under this formulation, the birth rate is constant , the death rate is linear , and the killing rate is linear This class of processes was studied by van Doorn and Zeifman (2005), who demonstrated a different approach for calculating . This alternative approach (not shown) confirms our derivation of (A4).

### Derivation of Asymptotic Bias

We are interested in estimating a single relative population size *η* (defined relative to the true population size, *N*), which must be incorporated into the transition density function at recombination sites under the SMC and SMC′. Under the SMC, this transition density function is (A9)This is equivalent to the conditional density (S6 in File S1) with the addition of a relative population size parameter. Under the SMC′, the transition function is (A10)which is equivalent to the conditional density (S3 in File S1) with a relative population size parameter included.

Under the SMC, given the local coalescence time *t*, the distance along the chromosome until the nearest recombination event (measured in units of *ρ*) is exponentially distributed with rate *t* (McVean and Cardin 2005). The likelihood function for a single relative population size *η* under the SMC is thus (A11)Under the SMC′, the likelihood function for a relative population size *η* is (A12)where is the exponential rate of encountering recombination events that change the coalescence time when the local coalescence time is *t* (Eriksson *et al.* 2009). Note that under the SMC, the length of a segment is independent of the relative population size *η* given the local coalescence time . This is not true for the SMC′, since the probability that the coalescence time changes at a recombination site depends on the population size.

As the length of the chromosome increases and the number of coalescence-time changes goes to infinity, the asymptotic maximum-likelihood estimate of the relative population size under the SMC is

(A13)Here the penultimate equality holds only if there is ergodic (*i.e.*, law-of-large-numbers-like) convergence of the sequence of pairs of coalescence times on either side of a recombination site under the ARG. In File S1, we show that the continuous-chromosome pairwise ARG is ergodic. That is, the mean coalescence time across a long chromosome converges to the mean coalescence time at a single point along the chromosome. We are unable to prove the ergodicity of the sequence of pairs of coalescence times at recombination sites where the coalescence time changes; instead, we note that (A13) is supported by simulation (see above). We also note that Wiuf (2006) proved the ergodicity of the discrete-locus ARG under a variety of neutral demographic models. A similarly in-depth proof may also apply for continuous-chromosome models, but we do not explore the point further.

In (A13), the ultimate equality follows from the fact that the joint distribution of coalescence times is marginally the same at recombination sites under the ARG and the SMC′. Numerical maximization of the double integral shows that the maximum-likelihood estimate of a single population size *N* under the pairwise SMC is asymptotically biased, with the asymptotic estimate being ∼.

Under the ARG, the stationary distribution of lengths between recombination events that change the local coalescence time (*i.e.*, the identity-by-descent segment length distribution) is slightly different from that of the SMC′. (They are different because subsequent recombination events “heal” with slightly different probabilities under the ARG, while under the SMC′, each subsequent recombination event heals with the same probability.) Under the ARG, the identity-by-descent (IBD) length distribution is not currently known. Given that under the SMC′ the maximum-likelihood estimator for a relative population size involves the observed lengths, it is not currently possible to calculate the asymptotic bias of the pairwise SMC′ maximum-likelihood estimator of a single population size. However, the IBD length distribution under the ARG is approximated very closely by the SMC′ IBD length distribution (Carmi *et al.* 2014), so the SMC′ estimator is likely to be nearly asymptotically unbiased.

## Footnotes

Supporting information is available online at http://www.genetics.org/lookup/suppl/doi:10.1534/genetics.114.173898/-/DC1.

*Communicating editor: R. Nielsen*

- Received December 17, 2014.
- Accepted March 12, 2015.

- Copyright © 2015 by the Genetics Society of America