simulating an admixture pulse in the recent past over a background of isolation-by-distance migration model

[VivaswatS]

Hello devs, thanks for the awesome software! I thought I'd post my methodology here and get feedback on the correctness of my approach and/or potential improvements. I'm using version 1.2.0

Sorry for the length of the post! The code to reproduce my simulations is presented at the end. Happy to provide more information.

I am currently working on a framework to model long-range migration events over a background of isolation-by-distance for a simple grid-like structure below (assuming equal population sizes in each deme & equal, symmetric migration rates between demes). Over this 8x10 grid, I add a long-range MassMigration event (pulse of individuals moving from source ($s$) to destination ($d$) forward in time with some strength $c$ and at some time in the recent past $\tau$, along the indicated arrow. I believe MassMigration is exactly the kind of behavior I want since: 1) I would like to keep the underlying migration rate map constant following the event, and 2) I want the source & destination demes to continue to exist as usual after the event.

Now, I want to be able to model the expected pairwise coalescent times post admixture as a function of the equilibrium pairwise coalescent time (null with no admixture) and the two parameters, $c$ and $\tau$.
To do this, I need to know: 1) does setting time=1 in add_mass_migration() actually simulate a migration event $1$ generation ago (and not $1/2N_e$ generations ago)?, 2) how does msprime simulate the movement of lineages between different demes from generation-to-generation given a migration matrix $M$, where $M_{ij}$ is the expected number of individuals scaled by population size moving from deme $i$ to deme $j$? I am modeling the probability that a lineage in deme $i$ ends up in deme $j$ in $\tau$ generations as $\exp((M-\text{diag}(M1))\tau)$ (maybe this is wrong?)

Currently, my model produces unbiased estimates for all expected pairwise coalescent times for $c \in [0,1]$ and $\tau \approx 0$ (in simulations, I see time $=0.00001$ to capture this instantaneous pulse), but overestimates $T_{dd}$ only for $\tau \in {1, 2, \ldots}$ (see figure below for estimated vs simulated values for within-deme coalescent times). This makes me wonder whether I'm simulating what I'm supposed to be modeling (maybe there's a change in units I'm missing) and whether, my model for movement of lineages is correct. But, there's a very slight chance, I could be calculating the pairwise coalescent times from the tree sequences wrong, so I post the code at the end to see if there are any glaring errors.

I use the following code (runs in about ~30 seconds on my machine):

import msprime
population_configurations = [
            msprime.PopulationConfiguration(sample_size=30, initial_size=1) for i in range(80)
        ]
   
# migration matrix (find code to construct matrix here: https://github.com/VivaswatS/feems/blob/main/docs/notebooks/adj_matrix.py)
#      [[0.   , 0.001, 0.   , ..., 0.   , 0.   , 0.   ],
#       [0.001, 0.   , 0.001, ..., 0.   , 0.   , 0.   ],
#       [0.   , 0.001, 0.   , ..., 0.   , 0.   , 0.   ],
#       ...,
#       [0.   , 0.   , 0.   , ..., 0.   , 0.001, 0.   ],
#       [0.   , 0.   , 0.   , ..., 0.001, 0.   , 0.001],
#       [0.   , 0.   , 0.   , ..., 0.   , 0.001, 0.   ]]
# ensuring Ne x m = 1 -> m = 0.001 is low making it both a rate & a probability (?)
demography = msprime.Demography.from_old_style(population_configurations=population_configurations,
                                                   migration_matrix=adj_matrix,
                                                   ignore_sample_size=True, 
                                                   Ne=1000)

c = 0.6; tau = 10; s = 23; d = 56

# individuals from s->d forward in time is lineages from d->s backward in time
demography.add_mass_migration(time=tau, source=d, dest=s, proportion=c)

ts = msprime.sim_ancestry(
        samples=[msprime.SampleSet(15,i,ploidy=2) for i in range(80)],
        demography=demography,
        sequence_length=1,
        recombination_rate=0,  # No recombination
        record_migrations=False, record_full_arg=False,
        model=[msprime.DiscreteTimeWrightFisher(duration=tau+1), msprime.StandardCoalescent()], # is this the right units? 
        num_replicates=1000
    )

# storing the tree sequences (not a bottleneck: doing it sequentially due to other code in for loop)
for i, tree_sequence in enumerate(ts):
            tree_sequence.dump("tree{:d}_8x10_c{:.1f}_t{:d}.tree".format(i,c,tau))

Then, I read the tree sequences back into python to construct the expected pairwise coalescent time matrix (takes about ~45 mins):

ndemes = 80; ninds = 2*15; nsnps = 1000

def get_coaltimemat(c, t):
    mytree = []
    # reading each tree into a list
    for j in np.arange(nsnps):
        mytree.append(tskit.load("tree{:d}_8x10_c{:.1f}_t{:d}.tree".format(j,t)))

    # creating a matrix to store the pairwise coalescent times (of size 80x80)
    coaltimemat = np.zeros((ndemes,ndemes))

    # looping through each SNP i.e. replicate
    for snp in range(nsnps):
        if snp % 100 == 0:
            print(f'Finished SNP{snp}...')

        tt = mytree[snp].first()
        # looping through each pair of individuals across all demes
        for i in range(ninds*ndemes):
            for j in range(i+1,ninds*ndemes):
                # storing the result in the appropriate deme index
                coaltimemat[int(i/ninds),int(j/ninds)] += tt.tmrca(i,j)

    # normalizing by the total number of comparisons 
    coaltimemat[np.diag_indices_from(coaltimemat)] /= np.sum(range(ninds))
    coaltimemat[np.triu_indices_from(coaltimemat,1)] /= ninds*ninds

    # symmetric matrix 
    coaltimemat = coaltimemat + coaltimemat.T
    coaltimemat[np.diag_indices_from(coaltimemat)] *= 0.5

    # normalizing by the total number of replicates
    coaltimemat /= nsnps

    plt.imshow(coaltimemat, cmap='Wistia'); plt.colorbar(); 
    return coaltimemat

coaltimemat = get_coaltimemat(c=0.6, t=10)

[petrelharp]

Hello, @VivaswatS! I've not gone into detail in your code, but I can answer your questions:

1. does setting time=1 in add_mass_migration() actually simulate a migration event 1 generation ago (and not 1/2Ne generations ago)?

It simulates the event 1 time unit ago, and we strongly suggest setting effective population sizes to be in actual numbers of individuals, so that everything is in real (i.e., unscaled) units; in which case one time unit is one generation. (However, I see you've set the initial_size of each deme to be 1; maybe you mean these to be of size 1?) In other words, msprime does not do any secret scaling under the hood; time units are generations in whatever size populations you've specified.

2. how does msprime simulate the movement of lineages between different demes from generation-to-generation given a migration matrix ?

Well, I'm hoping this question is answered here: https://tskit.dev/msprime/docs/latest/demography.html#sec-demography-migration - but, I think the expression you have there is correct.

I might also suggest simulating just one simulation with a longer sequence (and positive recombination rate), then just using ts.divergence(..., mode="branch"), which gives the mean pairwise coalescence time (multiplied by two), rather than simulating lots of separate trees and looping over samples to get the tmrcas. (It should be a lot faster, and have a lot less code.)

[jeromekelleher]

One minor point to add to @petrelharp's answer:

population_configurations = [
            msprime.PopulationConfiguration(sample_size=30, initial_size=1) for i in range(80)
        ]

I would recommend against mixing the deprecated old-style PopulationConfiguration API here and the newer Demongraphy approach, as it will potentially do confusing things. For example, I don't think the sample_size=30 has any effect here. You should be able to do the same thing here with demography = msprime.Demography.isolated_model([1] * 80), and then set your migration matrix afterwards.

[VivaswatS]

Thanks for that quick reply! It seems like I'm not overwriting the initial_size from the PopulationConfiguration when setting Ne in demography.from_old_style, but I assumed I was. If that's the case, I think that solves my issue!

I agree that I shouldn't be mixing the two, so I'll switch over to specifying the Demography object separately.