When migration leaves a clean trace: Decoupling migration from coalescence in the structured serial coalescent

When migration leaves a clean trace: Decoupling migration from coalescence in the structured serial coalescent | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results When migration leaves a clean trace: Decoupling migration from coalescence in the structured serial coalescent View ORCID Profile Hao Shen , John Novembre doi: https://doi.org/10.1101/2025.10.10.681523 Hao Shen 1 Department of Human Genetics, University of Chicago Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Hao Shen John Novembre 1 Department of Human Genetics, University of Chicago 2 Department of Ecology and Evolution, University of Chicago Find this author on Google Scholar Find this author on PubMed Search for this author on this site Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF Abstract With the rapid accumulation of population genomic data across space and time, there is an urgent need for demographic inference methods that incorporate explicit time-series modeling, achieve fine spatial resolution, and ensure clear identifiability between migration and coalescence rates. To address this need, we investigate pairwise genealogical processes under the structured serial coalescent, deriving evolution equations for pairwise branch length distributions and related statistics. By classifying the resulting relationships according to their parameter dependencies and computational complexity, we identify a class that is not only computationally tractable but also determined solely by migration rates. Building on this theoretical basis, we propose an inference framework for fine-resolution, time-varying migration rates inference and demonstrate its feasibility through simulation. We further outline how this framework can be extended to the joint estimation of migration and coalescence rates. Introduction Since Kingman introduced the coalescent ( Kingman, 1982a , b ), it has been widely used in population genetic theory and inferences. However, the standard coalescent theory typically oversimplifies reality by assuming contemporaneous sampling and panmixia, whereas empirical samples can come from different time periods and geographically structured populations with variable lineage migration and coalescence rates. These limitations have motivated key extensions: the serial coalescent for temporal stratification ( Rodrigo and Felsenstein, 1999 ) and the structured coalescent for population subdivision ( Notohara, 1990 ; Takahata, 1991 ). For models combining both temporal and spatial structure, terminology varies in the literature. Some authors retain “structured coalescent” ( Müller et al., 2017 ) while others use “structured serial coalescent” ( Ewing and Rodrigo, 2007 ). For clarity and distinction, we adopt the latter term throughout this paper. These extensions of the standard coalescent have enabled powerful demographic inference tools, particularly for estimating coalescent effective population sizes (e.g., inverse of the instantaneous local coalescence rates) and migration rates. For the serial coalescent, approaches using based on full gene genealogies using tools from Bayesian phylogentetics, were developed to estimate effective population sizes from temporally sampled data ( Drummond et al., 2002 , 2005 ). For the structured coalescent and structured serial coalescent, related methods using the full gene genealogy have been proposed to jointly infer migration rates and effective population sizes ( Beerli and Felsenstein, 2001 ; Ewing et al., 2004 ; Volz, 2012 ; Vaughan et al., 2014 ; De Maio et al., 2015 ; Müller et al., 2017 , 2018 ). However, because these methods rely on information of the full genealogy, their computational scalability with respect to the number of lineages and demes remains limited—even with modern algorithmic advances. Consequently, demographic inference for systems with hundreds or thousands of demes remains impractical under these frameworks. One effective way to solve the scalability problem is to use pairwise coalescence information instead of the whole genealogy. In the context of the structured coalescent, this means developing methods based on pairwise coalescence times. The distribution of pairwise coalescence times depends on both the migration rates and coalescence rates, while also being directly connected to empirical data through pairwise summary statistics. This makes pairwise coalescence times a useful basis for demographic inference. For example, the methods EEMS ( Petkova et al., 2016 ), FEEMS ( Marcus et al., 2021 ), and FRAME ( Shen and Novembre, 2025 ) all exploit the connection between expected pairwise coalescence times and the sample covariance matrix ( McVean, 2009 ). Meanwhile, the method MAPS ( Al-Asadi et al., 2019 ) uses the fact that the length distribution of long pairwise shared coalescent segments (LPSC segments, also know as the identity-by-descent or IBD tracts) is determined by the distribution of pairwise coalescence times ( Palamara et al., 2012 ; Carmi et al., 2013 ). All of these methods are scalable to hundreds of demes. While earlier methods like EEMS, FEEMS, and MAPS relied heavily on the assumption of symmetric migration, work by ( Lundgren and Ralph, 2019 ) and by ( Shen and Novembre, 2025 ) with the FRAME method now enable fine-resolution inference for cases with asymmetric gene flow. Despite these methodological advances, limitations persist. First, research across many domains requires time-series modeling to accommodate serial sampling and time-varying migration rates. In conservation genetics, samples collected at multiple time points can reveal how gene flow patterns change through time. In epidemiology, serially sampled viral genomes can be used to reconstruct transmission dynamics. In human evolutionary biology, ancient DNA provides genetic snapshots from different epochs, enabling reconstruction of migration histories within each period. FRAME, however, assumes migration–drift equilibrium under the structured coalescent and therefore cannot handle serial sampling or time-varying parameters. While it is possible to analyze data by dividing it into temporal slices and fitting each slice with an equilibrium model (as done demonstrated in the original MAPS and FRAME papers), it remains suboptimal compared to direct modeling under the structured serial coalescent framework. Second, these methods can have difficulty in distinguish the effects of migration rates and coalescence rates, leading to identifiability issues. (MAPS is an exception, but its effectiveness in asymmetric migration scenarios remains uncertain.) These limitations motivate the need for a new demographic inference method that can: (1) preserve the scalability advantages of pairwise coalescence-based approaches, (2) directly incorporate serial sampling and time-varying parameters under the structured serial coalescent framework, and (3) effectively decouple the inference of migration rates from that of coalescence rates. In this paper we demonstrate that such a method is theoretically achievable. Unlike the structured coalescent, where the equation for expected pairwise coalescence times in the structured coalescent was established ( Strobeck, 1987 ) even before the theory’s formal development, the analogous theory for pairwise branch lengths in the structured serial coalescent - defined as the summed branch lengths from two sampled lineages to their most recent common ancestor (the serial sampling counterpart to pairwise coalescence times) - remains largely undeveloped. To bridge this theoretical gap and facilitate new inference methods, we develop the foundational theory for pairwise branch length dynamics under the structured serial coalescent framework. Specifically, we derive evolution equations governing the probability density functions (PDFs) of pairwise branch lengths. By solving these equations, we uncover fundamental relationships among the PDFs and systematically classify them into distinct categories based on their parameter dependencies and computational complexity. This classification framework extends naturally to other quantities derived from pairwise branch lengths, as we demonstrate through parallel analysis of mean pairwise branch length dynamics and length distribution dynamics of LPSC segments. The general applicability of this classification enables evaluation of each relationship class’s potential for demographic inference. Crucially, we identify a promising class that exhibits exclusive dependence on lineage migration rates with no influence from coalescence parameters, and computational tractability scaling as O ( d 3 ) where d represents the number of demes. These properties directly address the fundamental challenges identified earlier - they resolve the migration-coalescence identifiability problem while maintaining the scalability required for fine-resolution analysis. This theoretical breakthrough establishes a foundation for two inference approaches. First, for studies focusing specifically on gene flow dynamics, these relationships enable highly scalable migration rate estimation that is decoupled from coalescence parameters. Second, joint estimation of migration and coalescence rates may also benefit from these relationships. A two-step approach—first estimating migration rates using these classes, then inferring coalescence rates conditioned on the migration parameters—could improve both identifiability and scalability, though the latter gain may be less pronounced than in migration-only inference. We outline inference frameworks for both approaches through a proof-of-concept example and demonstrate the feasibility of the first approach through preliminary simulation studies. The model The structured serial coalescent process models the combined effects of backward migration between demes and coalescence events occurring when lineages reside in the same deme. We represent backward migration as a continuous-time jump process on a weighted directed graph with nodes {1, 2, …, d }, each corresponding to a deme. The weight on the edge from node i to node j is denoted m ij ( t ), representing the backward migration rate from i to j at time t . Biologically, m ij ( t ) is the rate at which deme i receives ancestry from deme j . The transition rate matrix of this jump process is Q ( t ), and the associated Laplacian is defined as L ( t ) = − Q ( t ). The migration rates are encoded in the weighted adjacency matrix M ( t ) = Q ( t ) − diag Q ( t ). The coalescence rate in deme i is denoted by γ i ( t ), and all coalescence rates are collected into the vector γ ( t ) = ( γ i ( t )). In our model, we study the pairwise genealogical process using the pairwise branch length, defined as the sum of the branch lengths of the two lineages up to their MRCA. This measure generalizes pairwise coalescence time by explicitly accounting for temporal offsets between samples and can be viewed as the distance of the two lineages in terms of branch length. For an illustration of how pairwise branch length is computed, see Figure 1 . Download figure Open in new tab Fig. 1. Pairwise branch length. A toy example showing how pairwise branch lengths are computed. Assume we have three samples X, Y , and Z sampled at times t X , t Y , and t Z , respectively. The pairwise branch length between X and Y is . Similarly, the pairwise branch length between X and Z is , and between Y and Z is . Now let denote the random variable representing the pairwise branch length between two lineages, one sampled from deme i at time x and the other from deme j at time y , where x ≤ y and time is measured backward from the present ( t = 0). Let be the corresponding cumulative distribution function ; be its probability density function; be the matrix of the probability density functions (some times we will neglect b and just write f x,y to refer to the function). All the above functions have support in [ y − x , + ∞) and will be 0 for b < y − x . In the next three sections, we will explore the relationships among the probability density functions as well as two related statistics (expected pairwise branch length and survival function of LPSC segment length). Relationships among probability density functions of pairwise branch length We consider a toy example with three epochs: [ t 0 , t 1 ), [ t 1 , t 2 ), and [ t 2 , ∞), which we label as epoch 0, epoch 1, and epoch 2, respectively. The migration and coalescence rates are assumed to be piecewise constant in the first two epoches, that is, L ( x ) = L 0 , γ ( x ) = γ 0 for x ∈ [ t 0 , t 1 ) and L ( x ) = L 1 , γ ( x ) = γ 1 for x ∈ [ t 1 , t 2 ) (The migration and coalescence rates in the last epoch are not specified, as they are not used in this section). We now investigate relationships among the following time-ordered pairs of pairwise branch length pdf matrices: . By abuse of notation, we will also use these pairs to denote the corresponding relationships. These 8 relationships correspond to the 8 edges in Fig. 2 . Download figure Open in new tab Fig. 2. Graphical representation of the relationships among pairwise branch length probability density functions Graphically, these 8 edges in Fig. 2 fall naturally into three categories: the vertical edges, the diagonal edges and the horizontal edges. In what follows, we show that this graphical classification corresponds to a classification of the relationships themselves in terms of representational and computational complexity. We first examine the relationship , which is represented by a vertical edge. To this end, we ∈ study the evolution of for x ∈ [ t 0 , t 1 ). As we move from x to x + Δ t , coalescence does not occur, and only migration of the first lineage is relevant. This yields the following equation Expanding and ignoring higher-order terms, we cancel Δ t on both sides and obtain In matrix form, this partial differential equation becomes with the boundary condition of known . This PDE has the explicit solution Taking the limit as x → t 0 , we obtain the relationship . Following similar logic, the relationships and , which are likewise represented by vertical edges, can also be computed. These three relationships can be jointly expressed as These relationships are particularly valuable for two reasons. First, they depend only on backward migration rates and depend on L 0 (we represent the corresponding edges in blue in Fig. 2 ), while depends on L 1 . This makes them a promising tool for disentangling migration dynamics from coalescence effects (or effective population sizes) in complex migration–drift scenarios. Second, the computational complexity of these relationships is at most O ( d 3 ), and can often be reduced further by exploiting the (potentially) sparse structure of the Laplacian matrices. This scalability makes it possible to perform gene flow or migration inference at an unprecedented fine resolution, involving hundreds or even thousands of demes. As we shall see, the class of relationships represented by the vertical edges is the only one that possesses these favorable properties, and should therefore be fully exploited to facilitate inference. Now we turn to the relationships , which corresponds to a diagonal edge in Figure.2 . Let f x,b = f x,x,b denote the density matrix when the two lineages are sampled at the same time. For x ∈ [ t 0 , t 1 ), the associated partial differential equation is with boundary condition of known and f x ,0 = diag { γ 0 } for x ∈ [ t 0 , t 1 ). This equation is derived in a manner similar to equation (3) , with the key difference being that both lineages can migrate independently, and coalescence may occur when they are in the same deme. The PDE remains linear once we vectorize both sides. Let E i denote the matrix with a 1 in the ( i, i )-th entry and zeros elsewhere, and let ϵ i = vec( E i ) be its vectorization. Define , where ⊗ denotes the Kronecker product, define and . Then the solution to equation (6) is given by Taking the limit as x → t 0 , we obtain the relationship . A similar analysis applies to the relationship , which is represented by another diagonal edge. Together, the relationships represented by the diagonal edges can be written jointly as where . Before turning to the relationships represented by the horizontal edges, we would like to note two points. First, when the two lineages evolve together in f x,x,b , the model simplifies to a structured coalescent model. The probability density function of the pairwise coalescence times can be derived using the fact that the pairwise branch lengths are twice the pairwise coalescence times in the structured coalescent. Second, as t 2 goes to infinity in equation (8b, we obtain the stationary distribution of pairwise branch length under a structured coalescent model with migration rates coded by L 1 and coalescence rates coded by γ 1 . The stationary distribution of pairwise branch length distribution satisfies where . Because S 0 and S 1 depend on both migration and coalescence rates, the relationships represented by the diagonal edges are generally more complex and computationally intensive than those represented by the vertical edges. However, as we shall see next, they remain relatively straightforward to derive and interpret compared with the relationships represented by the horizontal edges, and thus serve as useful candidates for inferring coalescence rates or effective population sizes. For the relationships represented by the horizontal edges, namely , and , it is not easy to directly formulate and solve the corresponding PDEs. Since we have already established the relationships corresponding to the vertical and diagonal edges, we derive those for the horizontal edges indirectly (for details of the derivation, see Supplementary Information A), which yields These equations are highly complex both in terms of representation and computation. The difficulty arises from the fact that when we trace lineage 2 in f x,y,b along y , the process is not purely migrational: we must also account for the location of lineage 1 at the same time point. If lineage 1 happens to be in the same deme as lineage 2, there is a positive probability that they coalesce. Moreover, for lineage 1 to evolve from time x to time y , it must pass through all changes in the migration rates during the interval [ x, y ]. This is why the relationship , represented by the black horizontal edge in Fig. 2 , is particularly intricate, as it depends simultaneously on L 0 , L 1 , and γ 1 . Since all relationships represented by the horizontal edges can be effectively derived from those represented by the vertical and diagonal edges, and are typically much more complicated than the other two categories, it is important to avoid using them for inference. In a demographic inference framework based on pairwise branch lengths, one should rely only on relationships represented by vertical edges when estimating backward migration rates, and use both vertical and diagonal edges for the joint inference of migration rates and coalescence rates. Relationships among expected pairwise branch lengths The relationships among the probability density functions (pdfs) of pairwise branch lengths naturally extend to summary statistics derived from them. For example, consider the expected pairwise branch length . In the structured coalescent framework, the analogous quantity - expected pairwise coalescence times - are directly related to sample covariance structure ( McVean, 2009 ), a connection exploited by methods including EEMS ( Petkova et al., 2016 ), FEEMS ( Marcus et al., 2021 ), and FRAME ( Shen and Novembre, 2025 ). As we demonstrate in Supplementary Information D, such connection can extend to serially sampled data with slight modifications. Furthermore, advances in tree sequence reconstruction now allow for the direct estimation of from inferred branch lengths, though this requires careful consideration of the demographic assumptions inherent in the reconstruction process. Here, we establish the relationships among expected pairwise branch lengths in the structured serial coalescent. Similar to last section, we study a three-epoch toy model. The 8 relationships we would like to study are . These relationships corresponds to edges in Fig. 3 . Download figure Open in new tab Fig. 3. Graphical representation of the relationships among expected pairwise branch lengths We now study . When x ∈ [ t 0 , t 1 ), this expectation satisfies the partial differential equation with boundary condition of known . Using the identities L 1 d×d = 0 and for any α , the solution is given by Taking the limit as x → t 0 , we get the relationship . Following the same logic, the relationships represented by other vertical edges can also be computed. These relationships can be jointly written as The relationships corresponding to the diagonal edges can also be derived directly. Let db denote the matrix of expected pairwise branch lengths when both lineages are sampled at time x For x ∈ [ t 0 , t 1 ), this satisfies the ordinary differential equation Vectorizing and solving this equation gives where is the equilibrium expected pairwise branch length under the parameter set ( L 0 , γ 0 ) and satisfies the equation Now letting x → t 0 and applying the same logic to the epoch [ t 1 , t 2 ), we obtain the joint equation for the relationships represented by the diagonal edges where satisfies The relationships corresponding to the horizontal edges can also be derived indirectly (for a detailed derivation, see Supplementary Information B), and can be written as As expected, the relationships corresponding to the vertical edges depend only on migration rates and are relatively easy to compute, whereas those corresponding to the horizontal edges are the most complex in terms of both representation and computation. Relationships among survival functions of LPSC segment lengths In classical coalescent as well structured coalescent, it is well known that the length distribution of long pairwise shared coalescence segments (LPSC segments, or identity-by-descent tracts) is determined by the distribution of pairwise coalescence time ( Palamara et al., 2012 ; Carmi et al., 2013 ) and can thus be used for demographic inference. This is the idea underlying MAPS ( Al-Asadi et al., 2019 ), which infers the migration surfaces and coalescence rates based on the LPSC length distribution. In structured serial coalescent, following the logic of the previous two sections, we establish the relationships among the survival functions of the LPSC segment lengths. Let ρ x,y,µ be the matrix of survival functions such that it’s ( i, j ) th element is the probability that a LPSC segment between two randomly sampled lineages from deme i at time x and deme j at time y has length larger than µ . The probability of a random LPSC segment has length larger than µ given the pairwise branch length b and recombination rate r is simply e − rbu . And we obtain Here again we employ the three epoch toy model, and we would like to investigate the 8 relationships . Again, these 8 relationships correspond to 8 edges in Fig. 4 . Download figure Open in new tab Fig. 4. Graphical representation of the relationships among survival functions of LPSC segment lengths The easiest way to study the relationships represented by the vertical edges here to directly write out the integration out and use the relationships among the probability density functions, instead of solving the PDEs. Writing and using equation (5a) , (5b) and (5c) yields For the relationships represented by the diagonal edges, we can also derive analogous expressions. Applying equation (8a) , (8b) , with , we obtain Finally the relationships corresponding to the horizontal edges can also be derived indirectly, which are given by A step-by-step derivation of the equations presented in this section can be found in the Supplementary Information C. Once again, we observe that the relationships corresponding to the vertical edges are the most well-behaved, while those corresponding to the horizontal edges are the most complex in terms of both representation and computation. Inference framework Here we study how the above relationships can be used to establish inference frameworks. Similar to the previous sections, we set up a proof-of-concept example by assuming that samples are available at three time points, t 0 , t 1 , and t 2 , and that the goal is to infer the ‘effective’ demographic parameters in [ t 0 , t 1 ) and [ t 1 , t 2 ) using the relationships among the expected pairwise coalescence times, i.e., the matrices. Although the matrices are not directly observable, they are naturally connected to sample variance structures and reconstructed tree sequences, as discussed in previous sections. The protocols below do not detail the establishment of these connections but instead present a higher-level inference framework. For convenience, we refer to sample variance structures and reconstructed tree sequences collectively as ‘corresponding data’, and on this basis propose two demographic inference approaches: one specialized for inferring migration rates, and another designed for the joint estimation of migration and coalescence rates. Pure migration rates inference For pure migration rates inference in the proof-of concept example, we follow a simple three-step protocol: Estimate the matrix from the corresponding data. This estimation can be straightforward if genetic samples are available from every deme. However, for cases involving vacant demes (those without samples), imputation becomes necessary. One strategy is to use spatial imputation, which leverages geographical proximity to infer values for the vacant demes based on the directly estimated submatrix of that represents the demes with samples. Alternatively, a model-based imputation approach can be employed. This method operates under the assumption that the stochastic process governing lineages before time t 2 was driven by a set of equilibrium migration rates ( L ∞ ) and coalescence rates ( γ ∞ ). In other words, it assumes , where the equilibrium matrix is defined by the equation: Under this assumption, established equilibrium-based methods, such as FRAME ( Shen and Novembre, 2025 ), can be used to first infer the parameters L ∞ and γ ∞ , which in turn provides an estimate for the complete matrix. Estimate and L 1 based on corresponding data and estimated using equation (13c) . From equation (13c) , we know is a function of and L 1 , since in the first step we already have the inferred , we can write solely as a function of L 1 . Consequently, a connection between L 1 and the corresponding data is mediated through . This established relationship allows for the estimation of L 1 . Notably, if vacant demes prevent the direct estimation of from corresponding data, a model based imputation step can be performed simultaneously with the estimation of L 1 . Estimate , and L 0 based on corresponding data and estimated using equation (5b) . Here the logic is similar to that in step 2. One thing to be mentioned here is that equation (5a) can also be used to help estimate L 0 if can be properly estimated from corresponding data. In practice people can either use a more convenient one or use both of them. Joint inference of migration rates and coalescence rates While the joint inference of migration and coalescence rates involves additional steps, the procedure for estimating migration rates remains consistent with the gene-flow-only protocol. Therefore, for these established steps, we will simply repeat the instructions without further explanation. New procedures specific to joint inference will be explained in detail. The protocol in the proof-of-concept example is as follows: Estimate the matrix from the corresponding data. Estimate and L 1 based on corresponding data and estimated using equation (13c) . Estimate and γ 1 based on corresponding data, estimated and estimated L 1 using equation (8b) . The idea here again is to use as a mediater to connect γ 1 to corresponding data. An important thing to mention is that we directly use the estimate of L 1 from last step instead of estimating it jointly with L 1 . Estimate , L 0 based on corresponding data, estimated and estimated using equation (5a) and (5b). The underlying logic is analogous to step 3 of the gene-flow inference protocol. The key distinction is that is now already estimated from the previous step, which simplifies the direct application of equation (5a) . Again, in practice people can either use a more convenient equation or use both of them. Estimate and γ 0 based corresponding data and estimated using equation (8a) . This step follows the same logic as step 3. While the proof-of-concept example demonstrates the core framework, it can be generalized to cases with more than two time slices and other data types, such as sample covariance structures or LPSC segment data. Simulation Here, we test the feasibility of our first inference approach—inferring pure migration rates—through simulation studies. Similar to the previous analysis, we consider a model with three epochs, [ t 0 , t 1 ), [ t 1 , t 2 ), and [ t 2 , t 3 ). In each epoch, we impose a different migration topology and assign random effective population sizes (i.e., coalescence rates). For each combination of topologies, we sample 10 haplotype lineages per deme per epoch and use msprime (Baumdicker et al., 2022) to simulate 10, 000 tree sequences. Parameters are then inferred from the pairwise branch length data of these simulated datasets as follows: Estimate the matrix from the simulated dataset, which we denote by . For simplicity, in this step we directly use the sample mean of pairwise branch length, i.e. , since we have samples from all the demes. Estimate and L 1 based on corresponding data and estimated using equation (13c) . From the last step, we already have , and we first obtain a preliminary estimate of by taking the sample mean of the pairwise branch lengths, denoted . We then infer L 1 by solving an optimization problem with objective , where is the relative error matrix of equation (13c) computed from and Ψ is a smoothness penalty on edge weights, and λ is a hyperparameter selected by cross-validation (see Supplementary Information E for details). Once we obtain , we refine our estimate of using Estimate , and L 0 based on corresponding data and estimated using equation (5b) . The procedure follows the same logic as in last step: we first obtain a direct sample-mean estimate of , and then infer L 0 by solving an analogous optimization problem. Finally, we refine using the resulting . The results of the simulation study are shown in Fig. 5 . Figures 5(a–c) present the ground-truth topologies for each epoch. In Fig. 5a , the sequence of topologies is as follows: during epoch [ t 0 , t 1 ], we have a topology of large-scale directionally migrating lineages, where lineages from the left boundary migrate backward to the right boundary; during epoch [ t 1 , t 2 ], the topology is large-scale spatially converging lineages, where lineages converge backward to the center; and during epoch [ t 2 , +∞), the topology consists of a mixture of small-scale patterns. We denote these by topology 1, topology 2, and topology 3, respectively. Thus, Fig. 5a corresponds to the topology sequence 1 → 2 → 3. Figures 5b and 5c show two other settings obtained by rotating this sequence (2 → 3 → 1 for Fig. 5b and 3 → 1 → 2 for Fig. 5c ). The detailed setting and visualization of simulation parameters are provided in Supplementary Information F and Supplementary Fig. 1–3. Download figure Open in new tab Fig. 5. Simulation results. The first row ( a–c ) presents the ground truth topology sequences across three epochs. We define three distinct migration topologies: topology 1 is large-scale directionally migrating lineages, topology 2 is large-scale spatially converging lineages, and topology 3 is a mixutre of small scale patterns. Panel (a) shows the sequence 1 → 2 → 3, panel (b) the sequence 2 → 3 → 1, and panel (c) the sequence 3 → 1 → 2. The second row ( d–f ) shows the corresponding inferred migration patterns for each of the three settings. Figures 5(d–f) show the inferred migration patterns. In all three settings, our method largely recovers the underlying structures, including many of the small-scale patterns. The only exception is the spatially diverging lineages in topology 3, which are less clearly recovered. This might be due to fact that the center of the spatially diverging lineages is visited much less than the center of the spatially converging lineages. In previous studies such as EEMS ( Petkova et al., 2016 ), FEEMS ( Marcus et al., 2021 ), and FRAME ( Shen and Novembre, 2025 ), migration rates could only be inferred in a relative sense. Here since we directly use the data of pairwise branch lengths, we can estimate the absolute migration rates. We provide a more detailed discussion of this issue in the discussion section. Discussion This study establishes an analytical foundation for pairwise genealogical processes under the structured serial coalescent. By deriving and solving evolution equations for pairwise branch length distributions and their expectations, we characterize fundamental relationships governing the interplay of migration and coalescence across time intervals. Our systematic classification identifies distinct relationship categories— represented graphically as edge classes in Figs. 2 - 4 —and evaluates their inference utility through parametric dependencies and computational complexity. This analysis reveals that relationships represented by vertical edges depend exclusively on migration rates, which is also validated in our analysis of relationships among length distributions of LPSC segments. This decoupling of spatial dynamics from coalescence effects provides a powerful tool to resolve the identifiability problem due to the interaction between migration and coalescence and enables fine-resolution inference of migration rates dynamics. Leveraging this decoupling, we propose inference frameworks that fully exploits the migration-exclusive nature of relationships represented by vertical edges. Our approach employs a forward-in-time sequential inference strategy: beginning with inferring/imputing the focal quantity related to pairwise branch lengths and parameters in the oldest time point, we sequentially reconstruct migration rates (and coalescence rates) using relationships represented by vertical (and diagonal) edges. The architecture naturally accommodates time-stratified sampling and can be adapted to alternative summary statistics. Despite these advances, several challenges remain for future work. First, the optimization and inference procedure we employ is straightforward and simple. While this may suffice for the proof-of-concept example, a more thorough investigation of the objective, penalty, and cross-validation procedures will be necessary to develop formal inference methods based on different statistics derived from pairwise branch lengths. What’s more, for expected pairwise branch length based inference, estimation accuracy depends critically on reliable estimates of expected pairwise branch lengths. In empirical applications, pairwise branch lengths are estimated from inferred genealogical trees rather than the true underlying trees. The gap between these two sources requires careful evaluation, since most existing tree inference methods assume panmixia and do not incorporate a formal migration model. In addition, imputation procedures used to address demes with no samples may introduce further errors. In a sequential inference framework based on expected pairwise branch lengths, these errors can propagate, and similar propagation issues can arise in SNP-based inference and LPSC-based inference, although the sources of error can differ. Addressing this problem will be particularly important for inference spanning many epochs. One possible strategy to reduce error propagation is a block-wise sequential approach, in which parameters are inferred jointly across multiple consecutive epochs rather than strictly one epoch at a time. Second, our proof-of-concept example does not address the issue of long-range migration. This challenge has already been noted in previous studies ( Petkova et al., 2016 ; Marcus et al., 2021 ; Shen and Novembre, 2025 ). Strictly local migration renders the migration-rate matrix sparse, and such sparsity can be leveraged for computational efficiency ( Marcus et al., 2021 ) and parameter identifiability ( Lundgren and Ralph, 2019 ). However, a strictly local network may be insufficient, as real populations may also have experienced gene flow from long-range connections. The difficulty is that these long-range sources are not known a priori ; if they were, one could simply include them in the network. The FEEMS paper ( Marcus et al., 2021 ) proposed two strategies to address this issue: a “greedy” algorithm akin to TreeMix ( Pickrell and Pritchard, 2012 ) and a combination of the graphical lasso ( Friedman et al., 2008 ) with graph Laplacian smoothing ( Wang et al., 2016 ). The former has already been explored in FEEMSmix ( Shastry et al., 2025 ), whereas the latter remains to be investigated. Third, to further accelerate the inferences, it is important to adopt efficient algorithms to compute high-dimensional matrix exponentials—particularly for matrices with specific sparse structures. The case of pure migration rate estimation requires exponentiating a d × d matrix ( e Lt ), while coalescence rate estimation involves a d 2 × d 2 matrix ( e St ). While directly computing e Lt is still acceptable, directly computing e St can be computationally expensive for large d . In this case, computation and approximation strategies that maximally exploit the structure of S can be very helpful. As the inference challenges outlined earlier are properly addressed, this framework may support powerful applications across evolutionary biology, epidemiology, and conservation. In human evolution, it can augment reconstruction of demographic history from ancient DNA, clarifying migration patterns and their connections to cultural and environmental change. In epidemiological research, it could be aid reconstructions of pathogen transmission dynamics across space and time, facilitating the monitoring of disease spread. In conservation genetics, it may help reveal how gene flow among natural populations shifts through time and how these changes are shaped by environmental pressures such as habitat loss, fragmentation, and climate change. Together, these contributions may provide new avenues to analyze population histories where spatial and temporal dimensions interact—from deep evolutionary timescales to contemporary ecological monitoring. Code Availability The code to reproduce all results is available at https://github.com/ShenHaotv/When-migration-leaves-a-clean-trace . Acknowledgements We would like to thank Ahmed Selim, Sherif Negm and Egor Lappo for helpful discussions and feedback on the manuscript. The research was supported by funding from NIH NIGMS grant R35-GM149521 to John Novembre. Footnotes Minor revisions to improve clarity, consistency, and framing of the text. https://github.com/ShenHaotv/When-migration-leaves-a-clean-trace References ↵ Al-Asadi , H. , Petkova , D. , Stephens , M. , and Novembre , J. ( 2019 ). Estimating recent migration and population-size surfaces . PLOS Genetics , 15 ( 1 ): e1007908 . OpenUrl Baumdicker , F. et al. ( 2022 ). Efficient ancestry and mutation simulation with msprime 1.0 . Genetics , 220 : iyab229 . OpenUrl CrossRef PubMed ↵ Beerli , P. and Felsenstein , J. ( 2001 ). Maximum likelihood estimation of a migration matrix and effective population sizes in n subpopulations by using a coalescent approach . Proceedings of the National Academy of Sciences of the United States of America , 98 ( 8 ): 4563 – 4568 . OpenUrl Abstract / FREE Full Text ↵ Carmi , S. , Palamara , P. F. , Vacic , V. , Lencz , T. , Darvasi , A. , and Pe’er , I. ( 2013 ). The variance of identity-by-descent sharing in the wright-fisher model . Genetics , 193 ( 3 ): 911 – 928 . OpenUrl Abstract / FREE Full Text ↵ De Maio , N. , Wu , C.-H. , O’Reilly , K. M. , and Wilson , D. ( 2015 ). New routes to phylogeography: a Bayesian structured coalescent approximation . PLoS Genetics , 11 ( 8 ): e1005421 . OpenUrl PubMed ↵ Drummond , A. J. , Nicholls , G. K. , Rodrigo , A. G. , and Solomon , W. ( 2002 ). Estimating mutation parameters, population history and genealogy simultaneously from temporally spaced sequence data . Genetics , 161 ( 3 ): 1307 – 1320 . OpenUrl Abstract / FREE Full Text ↵ Drummond , A. J. , Rambaut , A. , Shapiro , B. , and Pybus , O. G. ( 2005 ). Bayesian coalescent inference of past population dynamics from molecular sequences . Molecular Biology and Evolution , 22 ( 5 ): 1185 – 1192 . OpenUrl CrossRef PubMed Web of Science ↵ Ewing , G. , Nicholls , G. , and Rodrigo , A. ( 2004 ). Using temporally spaced sequences to simultaneously estimate migration rates, mutation rate and population sizes in measurably evolving populations . Genetics , 168 ( 4 ): 2407 – 2420 . OpenUrl Abstract / FREE Full Text ↵ Ewing , G. and Rodrigo , A. ( 2007 ). Estimating population parameters using the structured serial coalescent with bayesian MCMC inference when some demes are hidden . Evolutionary Bioinformatics , 2 ( 2 ): 227 – 235 . OpenUrl PubMed ↵ Friedman , J. , Hastie , T. , and Tibshirani , R. ( 2008 ). Sparse inverse covariance estimation with the graphical lasso . Biostatistics , 9 ( 3 ): 432 – 441 . OpenUrl CrossRef PubMed Web of Science ↵ Kingman , J. F. C. ( 1982a ). The coalescent . Stochastic Processes and their Applications , 13 ( 3 ): 235 – 248 . OpenUrl CrossRef ↵ Kingman , J. F. C. ( 1982b ). On the genealogy of large populations . Journal of Applied Probability , 19 : 27 – 43 . OpenUrl CrossRef PubMed ↵ Lundgren , E. and Ralph , P. L. ( 2019 ). Are populations like a circuit? comparing isolation by resistance to a new coalescent-based method . Mol. Ecol. Resour ., 19 ( 6 ): 1388 – 1406 . OpenUrl CrossRef PubMed ↵ Marcus , J. , Ha , W. , Barber , R. F. , and Novembre , J. ( 2021 ). Fast and flexible estimation of effective migration surfaces . eLife , 10 : e61927 . OpenUrl CrossRef PubMed ↵ McVean , G. ( 2009 ). A genealogical interpretation of principal components analysis . PLoS Genetics , 5 : e1000686 . OpenUrl PubMed ↵ Müller , N. F. , Rasmussen , D. , and Stadler , T. ( 2018 ). MASCOT: Parameter and state inference under the marginal structured coalescent approximation . Bioinformatics , 34 ( 22 ): 3843 – 3848 . OpenUrl CrossRef PubMed ↵ Müller , N. F. , Rasmussen , D. A. , and Stadler , T. ( 2017 ). The structured coalescent and its approximations . Molecular Biology and Evolution , 34 ( 11 ): 2970 – 2981 . OpenUrl CrossRef PubMed ↵ Notohara , M. ( 1990 ). The coalescent and the genealogical process in geographically structured populations . J. Math. Biol ., 29 : 59 – 75 . OpenUrl CrossRef PubMed Web of Science ↵ Palamara , P. F. , Lencz , T. , Darvasi , A. , and Pe’er , I. ( 2012 ). Length distributions of identity by descent reveal fine-scale demographic history . The American Journal of Human Genetics , 91 ( 5 ): 809 – 822 . OpenUrl CrossRef PubMed ↵ Petkova , D. , Novembre , J. , and Stephens , M. ( 2016 ). Visualizing spatial population structure with estimated effective migration surfaces . Nature Genetics , 48 : 94 – 100 . OpenUrl CrossRef PubMed ↵ Pickrell , J. K. and Pritchard , J. K. ( 2012 ). Inference of population splits and mixtures from genome-wide allele frequency data . PLoS Genetics , 8 ( 11 ): e1002967 . OpenUrl ↵ Crandall , K. A. , editor Rodrigo , A. G. and Felsenstein , J. ( 1999 ). Coalescent approaches to hiv population genetics . In Crandall , K. A. , editor, The Evolution of HIV , pages 233 – 272 . Johns Hopkins University Press . ↵ Shastry , V. , Musiani , M. , and Novembre , J. ( 2025 ). Jointly representing long-range genetic similarity and spatially heterogeneous isolation-by-distance . bioRxiv . https://doi.org/10.1101/2025.02.10.637386 . ↵ Shen , H. and Novembre , J. ( 2025 ). Fine-resolution asymmetric migration estimation . bioRxiv . https://doi.org/10.1101/2025.05.29.656894 . ↵ Strobeck , C. ( 1987 ). Average number of nucleotide differences in a sample from a single subpopulation: a test for population subdivision . Genetics , 117 : 149 – 153 . OpenUrl Abstract / FREE Full Text ↵ Takahata , N. ( 1991 ). Genealogy of neutral genes and spreading of selected mutations in a geographically structured population . Genetics , 29 : 585 – 595 . OpenUrl ↵ Vaughan , T. G. , Kühnert , D. , Popinga , A. , Welch , D. , and Drummond , A. J. ( 2014 ). Efficient Bayesian inference under the structured coalescent . Bioinformatics , 30 ( 16 ): 2272 – 2279 . OpenUrl CrossRef PubMed ↵ Volz , E. M. ( 2012 ). Complex population dynamics and the coalescent under neutrality . Genetics , 190 ( 1 ): 187 – 201 . OpenUrl Abstract / FREE Full Text ↵ Wang , Y.-X. , Sharpnack , J. , Smola , A. J. , and Tibshirani , R. J. ( 2016 ). Trend filtering on graphs . Journal of Machine Learning Research , 17 ( 105 ). View the discussion thread. Back to top Previous Next Posted October 11, 2025. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following When migration leaves a clean trace: Decoupling migration from coalescence in the structured serial coalescent Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share When migration leaves a clean trace: Decoupling migration from coalescence in the structured serial coalescent Hao Shen , John Novembre bioRxiv 2025.10.10.681523; doi: https://doi.org/10.1101/2025.10.10.681523 Share This Article: Copy Citation Tools When migration leaves a clean trace: Decoupling migration from coalescence in the structured serial coalescent Hao Shen , John Novembre bioRxiv 2025.10.10.681523; doi: https://doi.org/10.1101/2025.10.10.681523 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Evolutionary Biology Subject Areas All Articles Animal Behavior and Cognition (7618) Biochemistry (17636) Bioengineering (13860) Bioinformatics (41847) Biophysics (21401) Cancer Biology (18536) Cell Biology (25424) Clinical Trials (138) Developmental Biology (13353) Ecology (19860) Epidemiology (2067) Evolutionary Biology (24287) Genetics (15583) Genomics (22463) Immunology (17701) Microbiology (40300) Molecular Biology (17141) Neuroscience (88434) Paleontology (666) Pathology (2825) Pharmacology and Toxicology (4813) Physiology (7633) Plant Biology (15107) Scientific Communication and Education (2042) Synthetic Biology (4285) Systems Biology (9808) Zoology (2268)