Configuring demography

The Coalescent expects a demography object to be passed to it, which can be configured in various ways. When constructing a demography object, you can directly specify the time points at which the population sizes or migration rates change.

import phasegen as pg

d = pg.Demography(
    pop_sizes={'pop_0': {0: 1}, 'pop_1': {0: 2.5, 1: 0.8}},
    migration_rates={
        ('pop_0', 'pop_1'): {0: 1.7, 0.7: 2},
        ('pop_1', 'pop_0'): {0: 3}
    }
)

d.plot();

Alternatively, you can configure the demography object after construction by adding demographic events. Below, we add PopSizeChange and MigrationRateChange events to the demography object.

d = pg.Demography()

d.add_event(pg.PopSizeChange(pop='pop_0', time=0, size=1))
d.add_event(pg.PopSizeChange(pop='pop_1', time=0, size=2.5))
d.add_event(pg.PopSizeChange(pop='pop_1', time=1, size=0.8))

d.add_event(pg.MigrationRateChange(source='pop_0', dest='pop_1', time=0, rate=1.7))
d.add_event(pg.MigrationRateChange(source='pop_0', dest='pop_1', time=0.7, rate=2))
d.add_event(pg.MigrationRateChange(source='pop_1', dest='pop_0', time=0, rate=3))

d.plot();

This is similar to the Msprime demography API, and we can easily convert to an msprime.Demography object. Note that the reverse, converting an msprime demography to a native Demography object, is not currently supported due to phasegen’s inherent restriction to discrete rate changes.

d_msprime = d.to_msprime()

Discretizing continuous demographies

There are also utilities for discretizing continuous demographies. In the example below, we create a discretized demography by passing a continuous callback function to DiscretizedRateChange. You can freely combine this with other demographic events. Note that the total runtime of the computations is linear in the number of epochs, i.e., it is roughly a multiple of the number of epochs.

d = pg.Demography(
    events=[
        pg.DiscretizedRateChange(
            trajectory=lambda t: 1.5 + 0.5 * t,
            pop='pop_0',
            start_time=0,
            end_time=5,
            step_size=0.5
        )
    ]
)

d.plot();

For exponential growth or decline, you can also make use of ExponentialPopSizeChanges.

d = pg.Demography(
    events=[
        pg.ExponentialPopSizeChanges(
            initial_size={'pop_0': 1.5},
            growth_rate=0.5,
            start_time=0,
            end_time=8,
            step_size=0.5
        )
    ]
)

d.plot();

Population splits

Population splits (forwards in time) correspond to population mergers (backwards in time). Since phasegen does not support deterministic lineage movements due to its inherent structure, we can model a population split by specifying a large unidirectional migration rate from the derived to the ancestral population. Below, we model a population split where pop_0 splits from pop_1 at time 2. This corresponds to a population merger of the two populations at time 2 backwards in time, and we thus need to initialize both populations at time 0 in the present.

d = pg.Demography(
    pop_sizes={'pop_0': 1, 'pop_1': 3},
    events=[
        pg.PopulationSplit(
            derived='pop_0',
            ancestral='pop_1',
            time=2
        )
    ]
)

Plotting the migration rates, we see that there is a large migration rate from pop_0 to pop_1 at time 2. All migration rates to the derived population (pop_0) are furthermore set to 0 at the time of the split.

d.plot_migration();

Wrapping in a Coalescent object, by specifying the initial numbers of lineages in each population, we can visualize the tree height distribution. We see that the probability of absorption is 0 after to the split (forwards in time). This is because our scenario represents a clean population split, where the derived population is completely isolated from the ancestral population after the split, which makes coalescence between the two populations impossible.

coal = pg.Coalescent(
    n={'pop_0': 4, 'pop_1': 4},
    demography=d
)

coal.tree_height.plot_pdf();

Population mergers

Population mergers (forwards in time) correspond to population splits (backwards in time). Mergers of populations that were completely isolated prior to the merger are difficult to model in a coalescent framework. This is because, backwards in time, we would end up with isolated populations which would not be able to coalesce. We can thus only model mergers of populations, provided they are in contact again eventually. Below an example

import matplotlib.pyplot as plt

coal = pg.Coalescent(
    n={'pop_0': 8, 'pop_1': 0},
    demography=pg.Demography(
        pop_sizes={'pop_0': 1, 'pop_1': 1},
        events=[
            pg.MigrationRateChanges(
                {
                    ('pop_1', 'pop_0'): {2: 0, 8: 0.1},
                    ('pop_0', 'pop_1'): {0.5: 1, 3: 0, 8: 0.1}
                }
            )
        ]
    )
)

_, axs = plt.subplots(1, 2, figsize=(10, 3.5))
coal.demography.plot_migration(ax=axs[0])
coal.tree_height.plot_pdf(ax=axs[1]);

Note

Migration rates between demes are assumed to be 0 if not specified otherwise.