When is it appropriate to use add_mass_migration()?

[David-Peede]

Hello again!

I am interested in studying the coalescent properties of introgression, which for me involves comparing simulated results vs analytical solutions. However, I just noticed a small quirk that I didn't initially think about. In the legacy version of msprime I modeled introgression and population divergences using the MassMigration() demography object. For example:

import msprime

def one_pulse(
        seq_len,
        Ne_P1,
        Ne_P2,
        Ne_P3,
        Tp3,
        Tp2,
        Tgf,
        mu,
        gamma,
        r=None,
        num_seq=1,
        random_seed=None,
):
    """
    Simulates three populations ((P1, P2), P3) where the only
    demographic events are population divergences and one pulse of
    introgression. Sequence length is in base pairs and both mutation and
    recombination rate are in base pairs per generation. This function returns
    a tree sequence file.
    """
    # Define simulation parameters
    samp_size = num_seq # Number of sequences
    prob_introgress = gamma # Probability of a lineage introgressing
    # Define samples
    samples = [msprime.Sample(population=0,time=0)]*samp_size # P1
    samples.extend([msprime.Sample(population=1,time=0)]*samp_size) # P2
    samples.extend([msprime.Sample(population=2,time=0)]*samp_size) # P3
    # Configure populations with a constant population size
    population_configurations = [
        msprime.PopulationConfiguration(initial_size=Ne_P1),
        msprime.PopulationConfiguration(initial_size=Ne_P2),
        msprime.PopulationConfiguration(initial_size=Ne_P3),
        ]
    # Define demographic events
    demographic_events = [
        msprime.MassMigration(
                time=Tgf, source=1, destination=2, proportion=prob_introgress, # Introgression event
                ),
        msprime.MassMigration(
                time=Tp2, source=1, destination=0, proportion=1.0, # P1 and P2 merge
                ),
        msprime.MassMigration(
                time=Tp3, source=2, destination=0, proportion=1.0, # P12 and P3 merge
                ),
        ]
    # Simulate in msprime
    ts = msprime.simulate(
        samples=samples,
        length=seq_len,
        recombination_rate=r,
        population_configurations=population_configurations,
        demographic_events=demographic_events,
        random_seed=random_seed,
        )
    return ts

When "migrating" over to the new msprime I used the add_mass_migration() object to simulate the introgression event and then used the add_population_split() object to simulate the population divergences. For example:

import msprime

def IUA_pop_split(
    N,
    f,
    Tgf,
    Tp2,
    Tp3,
    donor='P3',
    recipient='P2',
    r=None,
    reps=None,
    sequence_length=None,
    random_seed=None,
):
    """
    Simulates three lineages ((P1, P2), P3) where the only
    demographic events are population divergences and one pulse
    of introgression. This function returns a tree sequence file.
    """
    # Intialize the demographic model.
    demography = msprime.Demography()
    # We assume constant and equal effective population sizes for
    # all lineages.
    demography.add_population(name='P1', initial_size=N)
    demography.add_population(name='P2', initial_size=N)
    demography.add_population(name='P3', initial_size=N)
    demography.add_population(name='P12', initial_size=N)
    demography.add_population(name='P123', initial_size=N)
    # Introgression from P3 to P2 with a probaility f at time Tgf.
    demography.add_mass_migration(
        time=Tgf, source=recipient, dest=donor, proportion=f,
    )
    # P1 and P2 merge into P12
    demography.add_population_split(
        time=Tp2, derived=['P1', 'P2'], ancestral='P12'
    )
    # P12 and P3 merge into P123
    demography.add_population_split(
        time=Tp3, derived=['P12', 'P3'], ancestral='P123'
    )
    # Simulate!
    ts = msprime.sim_ancestry(
        samples=[
            msprime.SampleSet(1, ploidy=1, population='P1'),
            msprime.SampleSet(1, ploidy=1, population='P2'),
            msprime.SampleSet(1, ploidy=1, population='P3'),
        ],
        recombination_rate=r,
        demography=demography,
        random_seed=random_seed,
        num_replicates=reps,
        
    )
    # As a sanity check let's also print the demograpghy debugger.
    #print(demography.debug())
    return ts

Then I came across @jeromekelleher and Konrad Lohse's introgression example in the tskit tutorials and noticed that they opted to use the add_mass_migration() to simulate the introgression event AND population divergences, which made me wonder "is there a reason why they used add_mass_migration() vs add_population_split() to simulate models of introgression"? So I then decided to write up another version of the IUA introgression model only using the add_mass_migration() object :

import msprime

# Define an IUA model of introgression.
def IUA_migration(
    N,
    f,
    Tgf,
    Tp2,
    Tp3,
    donor='P3',
    recipient='P2',
    r=None,
    reps=None,
    sequence_length=None,
    random_seed=None,
):
    """
    Simulates three lineages ((P1, P2), P3) where the only
    demographic events are population divergences and one pulse
    of introgression. This function returns a tree sequence file.
    """
    # Intialize the demographic model.
    demography = msprime.Demography()
    # We assume constant and equal effective population sizes for
    # all lineages.
    demography.add_population(name='P1', initial_size=N)
    demography.add_population(name='P2', initial_size=N)
    demography.add_population(name='P3', initial_size=N)
    # Introgression from P3 to P2 with a probaility f at time Tgf.
    demography.add_mass_migration(
        time=Tgf, source=recipient, dest=donor, proportion=f,
    )
    # P1 and P2 merge into P12
    demography.add_mass_migration(
        time=Tp2, source='P2', dest='P1', proportion=1,
    )
    # P12 and P3 merge into P123
    demography.add_mass_migration(
        time=Tp3, source='P3', dest='P1', proportion=1,
    )
    # Simulate!
    ts = msprime.sim_ancestry(
        samples=[
            msprime.SampleSet(1, ploidy=1, population='P1'),
            msprime.SampleSet(1, ploidy=1, population='P2'),
            msprime.SampleSet(1, ploidy=1, population='P3'),
        ],
        recombination_rate=r,
        demography=demography,
        random_seed=random_seed,
        num_replicates=reps,
        
    )
    # As a sanity check let's also print the demograpghy debugger.
    #print(demography.debug())
    return ts

So I then planned to compare the two models by calculating the average internal branch length out of 1,000,000 replicates that would generate an ABBA or BABA site and compare the two models against the analytical solutions provided by Durand et al, 2011. This process required a bit more code, so for brevity, I will just show the results but I am happy to share the code and output files if need be!

Tables: f represents the introgression proportion, E[X] represents the analytical solutions, Mig (X) represents simulations using the IUA_migration() function, and Split (X) represents simulations using the IUA_pop_split() function.

f	E[ABBA]	Mig (ABBA)	Split (ABBA)	E[BABA]	Mig (BABA)	Split (BABA)
0	365874	364756	365692	365874	365311	363459
0.25	879884	880202	880976	329884	331127	329866
0.5	1.39389e+06	1.39185e+06	1.39617e+06	293894	294847	294089
0.75	1.9079e+06	1.90795e+06	1.91186e+06	257904	257892	259057
1	2.42191e+06	2.42322e+06	2.4235e+06	221914	220563	224111

f	E[ABBA] - Mig (ABBA)	E[ABBA] - Split (ABBA)	E[BABA] - Mig (BABA)	E[BABA] - Split (BABA)
0	1118.56	181.983	563.74	2415.13
0.25	-317.971	-1091.8	-1242.53	18.7147
0.5	2039.9	-2276.85	-952.427	-195.076
0.75	-41.3862	-3953.87	12.3273	-1152.73
1	-1306.32	-1582.06	1351.33	-2196.77

Lastly, I made a plot to see if there were any stark differences.

So this is a long-winded way to ask when studying introgression does it make more sense to model it using the add_mass_migration()? Based on this very small example I would naively think yes, but I don't understand why this is or is not the case. And more generally what are examples of advanced features where add_mass_migration() is appropriate?

Thanks in advance and I apologize in advance if anything I said was unclear.
-Dave

[David-Peede]

Update: After writing my first post I realized I should probably also benchmark the model written in the legacy version of msprime and recreate the tables and figure including those simulations. The only difference here is now there is an extra category where Old (X) represents simulations using the one_pulse() function defined above.

f	E[ABBA]	Old (ABBA)	Mig (ABBA)	Split (ABBA)	E[BABA]	Old (BABA)	Mig (BABA)	Split (BABA)
0	365874	364756	364756	365692	365874	365311	365311	363459
0.25	879884	880202	880202	880976	329884	331127	331127	329866
0.5	1.39389e+06	1.39185e+06	1.39185e+06	1.39617e+06	293894	294847	294847	294089
0.75	1.9079e+06	1.90795e+06	1.90795e+06	1.91186e+06	257904	257892	257892	259057
1	2.42191e+06	2.42322e+06	2.42322e+06	2.4235e+06	221914	220563	220563	224111

f	E[ABBA] - Old (ABBA)	E[ABBA] - Mig (ABBA)	E[ABBA] - Split (ABBA)	E[BABA] - Old (BABA)	E[BABA] - Mig (BABA)	E[BABA] - Split (BABA)
0	1118.56	1118.56	181.983	563.74	563.74	2415.13
0.25	-317.971	-317.971	-1091.8	-1242.53	-1242.53	18.7147
0.5	2039.9	2039.9	-2276.85	-952.427	-952.427	-195.076
0.75	-41.3862	-41.3862	-3953.87	12.3273	12.3273	-1152.73
1	-1306.32	-1306.32	-1582.06	1351.33	1351.33	-2196.77

Since the model using add_mass_migration() is identical to the legacy version I am slightly more confident that is how I should model introgression, but I still don't understand why...

Thanks again!
-Dave

[jeromekelleher]

The "mass migration" is a power-user feature @David-Peede, and can be used to imitate the effects population split and admixture events. It's appropriate to use mass migrations whenever you're modelling something that isn't a population split or admixture, so the introgression situation is one of these occasions. I'm guessing the tskit tutorial doesn't use population splits because it was ported from the old code as simply as possible, and we didn't think of using the population split (PR welcome if you'd like to update it!)

As an aside, have you considered using Demes to construct your model?There is excellent tooling available for visualising models in demes and it's quite a bit more readable than the msprime notation.

[David-Peede]

Hi @jeromekelleher that makes sense and unfortunately I have yet to start "migrating" over to Demes but I will make that my next priority! Sorry again for the long thread about a relatively straightforward answer, but I really appreciate it!