Soft selection with discrete demes

This page describes implementing models of demographic events affecting multiple demes. This functionality was first released in version 0.6.0 and makes use of low-level types added in 0.5.3.

Note

As of 0.6.0, these features only apply to simulations using tree sequence recording.

Overview

This chapter describes how to generate detailed demographic models of individuals evolving withing discrete sub-populations, or demes.

If you want to skip the details and see basic models provided by fwdpy11, then see The demography debugger.

Once you’ve digested this section, you may want to read Different effect sizes of mutations in different demes.

The model

The model here is soft selection [Levene1953], meaning that the number of breeding individuals (“adults”) in each deme is fixed at a certain value. A nice overview of this model and how it compares to others in [Felsenstein1976]. You may also find [Christiansen1974] and [Christiansen1975] useful.

Each generation, we generate offspring (“juveniles”) in each deme. In the absence of migration, all parents of all offspring come from the offspring deme. With a migration matrix, we first choose the “source” deme, pick parents from that deme, and then create the offspring in the offspring deme. Thus, we are modeling juvenile migration. Selfing is a property of demes and is applied to parents: we choose a parental deme, then choose a parent, and then decide if that parent selfs.

Note

The migration behavior changed in 0.6.2! (This is mainly a note for the developers.)

The timings of events

Below, we discuss various events that may happen. These event types include things like deme size changes, “mass migration” events, etc.. These events will occur at a certain time in a simulation. That time refers to the birth time of a generation and the events are applied prior to generating offspring, meaning that the events happen to the parents. For example, if half of deme zero moves and colonizes a new deme (deme 1), then that means that half of the current alive individuals (possible parents) have their deme field changed from zero to one prior to generating any offspring.

The objects events all have an attribute called when, which parameterizes when the event occurs in a simulation. The value of when is with respect to the current generation time of the population (fwdpy11.DiploidPopulation.generation). When events are registered to occur prior to this time, you will see warnings. See fwdpy11.evolvets() for details.

Events scheduled for prior than the population’s current time are allowed because there are valid modeling reasons to allow them. For example, you may want to evolve for a while and then change some other model parameter like the recombination rate, and then keep evolving.

Setting the initial demes in a simulation

At the start of a simulation, you may assign diploids to demes when constructing an instance of fwdpy11.DiploidPopulation. For example, to initialize a population with 25 individuals in demes 0 and 1:

In [1]: pop = fwdpy11.DiploidPopulation([25, 25], 1.0)

In [2]: md = np.array(pop.diploid_metadata, copy=False)

In [3]: pop.deme_sizes()
Out[3]: (array([0, 1], dtype=int32), array([25, 25]))

In [4]: for m in pop.diploid_metadata:
   ...:     for n in m.nodes:
   ...:         assert m.deme == pop.tables.nodes[n].deme
   ...: 

Another method involves mass migration events at the beginning of a simulation. See Mass migrations.

The DiscreteDemography class

The demographic events are stored in instances of fwdpy11.DiscreteDemography. These events, whose interface is described below, are passed in list objects when creating a fwdpy11.DiscreteDemography instance.

These instances may be used to parameterize the demography field of a fwdpy11.ModelParams instance. To illustrate this, here is a function that we’ll use repeatedly below:

In [5]: def setup_and_run_model(pop, ddemog, simlen, recorder=None, seed=654321):
   ...:     pdict = {
   ...:         "nregions": [],
   ...:         "sregions": [],
   ...:         "recregions": [],
   ...:         "rates": (0, 0, 0,),
   ...:         "gvalue": fwdpy11.Multiplicative(2.0),
   ...:         "demography": ddemog,
   ...:         "simlen": simlen,
   ...:     }
   ...:     params = fwdpy11.ModelParams(**pdict)
   ...:     rng = fwdpy11.GSLrng(seed)
   ...:     fwdpy11.evolvets(rng, pop, params, 100, recorder)
   ...: 

We will also define a simple class to record all deme sizes over time:

# fmt: off
In [6]: class SizeTracker(object):
   ...:     def __init__(self):
   ...:         self.data = []
   ...:     def __call__(self, pop, sampler):
   ...:         self.data.append((pop.generation, pop.N, pop.deme_sizes()))
   ...: 

# fmt: on

Compatibility with previous versions of fwdpy11

Changed in version 0.8.0: You now must specify simlen manually.

Previous versions only supported size changes within a single deme. These size changes were parameterized via a numpy array specifying the size at each time point. It is still possible to specify the demography using that approach:

In [7]: N = np.array([10] * 10 + [5] * 5 + [10] * 10, dtype=np.uint32)

In [8]: pdict = {
   ...:     "nregions": [],
   ...:     "sregions": [],
   ...:     "recregions": [],
   ...:     "rates": (0, 0, 0,),
   ...:     "gvalue": fwdpy11.Multiplicative(2.0),
   ...:     "demography": fwdpy11.DiscreteDemography(set_deme_sizes=N),
   ...:     "simlen": len(N),
   ...: }
   ...: 

In [9]: params = fwdpy11.ModelParams(**pdict)

In [10]: rng = fwdpy11.GSLrng(654321)

In [11]: pop = fwdpy11.DiploidPopulation(10, 1.0)

In [12]: fwdpy11.evolvets(rng, pop, params, 100)

Event types

The following sub-sections describe the various types of demographic events allowed during a simulation.

Mass migrations

Mass migration events represent the “bulk” movement of individuals in a single generation. Such events allow you to model population splits, merges, etc..

These events are represented by instances of fwdpy11.MassMigration. Currently, you create instances of this type using one of the following two functions:

As the name implies, the first function creates an event that copies individuals from a source deme to a destination. The latter moves them.

Both functions take five arguments, which may be used either named or unnamed. In order, they are:

  • when: the time (generation) when the event will occur

  • source: the ID of the source deme

  • destination: the ID of the destination deme

  • fraction: the fraction (proportion) of source moved/copied to dest.

  • resets_growth_rate: If True, the event resets the growth rate to fwdpy11.NOGROWTH in both source and dest. If False, growth rates remain unchanged. The default is False.

Note

When a mass migration event copies individuals from deme, the individuals copied are sampled without replacement. Thus, if the fraction copied is 1.0, then every individual is copied.

These operations act on proportions of populations rather than on numbers of individuals. Multiple events in a single generation are allowed, see Multiple mass migrations.

Setting the initial state of a simulation

Let’s look at an example where we use mass migration events to set up “who is where” at the start of a simulation. Since events happen in the parental generation, we can use mass migrations to set up what demes individuals are in by applying events at generation 0.

The main difference between this method and that shown in Setting the initial demes in a simulation is that these events move or copy random individuals to new demes whereas using the __init__ approach builds the individuals in each deme sequentially.

For example, if we wish to start a simulation with 50 individuals in demes 0 and 50 in deme 1, we have two options:

  1. Start with 50 individuals and copy them to deme 1 in generation 0

  2. Start with 100 individuals and move half of them to deme 1 in generation 0

Here is the version implemented via a copy:

In [13]: pop = fwdpy11.DiploidPopulation(50, 1.0)

In [14]: copy = [fwdpy11.copy_individuals(when=0, source=0, destination=1, fraction=1.0)]

In [15]: ddemog = fwdpy11.DiscreteDemography(mass_migrations=copy)

In [16]: setup_and_run_model(pop, ddemog, 1)

In [17]: pop.deme_sizes()
Out[17]: (array([0, 1], dtype=int32), array([50, 50]))

Here is what our object looks like:

In [18]: copy[0]
Out[18]: fwdpy11.MassMigration(when=0, source=0, destination=1, fraction=1.0, move_individuals=False, resets_growth_rate=True)

Here is the version using a move:

In [19]: pop = fwdpy11.DiploidPopulation(100, 1.0)

In [20]: move = [fwdpy11.move_individuals(0, 0, 1, 0.5)]

In [21]: ddemog = fwdpy11.DiscreteDemography(mass_migrations=move)

In [22]: setup_and_run_model(pop, ddemog, 1)

In [23]: pop.deme_sizes()
Out[23]: (array([0, 1], dtype=int32), array([50, 50]))

For comparison, here is the object specifying the move:

In [24]: move[0]
Out[24]: fwdpy11.MassMigration(when=0, source=0, destination=1, fraction=0.5, move_individuals=True, resets_growth_rate=True)

Multiple mass migrations

To specify multiple events, simply add more events to your list. The events to not have to be sorted in any specific way. Any sorting requirements get handled internally.

Multiple events involving the same source population in the same generation need some explaining. If the events are copies, things will tend to “just work”:

In [25]: pop = fwdpy11.DiploidPopulation(50, 1.0)

In [26]: copy = [fwdpy11.copy_individuals(0, 0, 1, 1.0), fwdpy11.copy_individuals(0, 0, 2, 1.0)]

In [27]: ddemog = fwdpy11.DiscreteDemography(mass_migrations=copy)

In [28]: setup_and_run_model(pop, ddemog, 1)

In [29]: pop.deme_sizes()
Out[29]: (array([0, 1, 2], dtype=int32), array([50, 50, 50]))

When the events are moves, it is not possible to move more than 100% of the individuals. Attempting to do so will raise a ValueError exception:

In [30]: pop = fwdpy11.DiploidPopulation(50, 1.0)

# Move all of deme 0 into demes 1 and 2,
# which means we're trying to move 200%
# of deme 0...
In [31]: move = [fwdpy11.move_individuals(0, 0, 1, 1.0), fwdpy11.move_individuals(0, 0, 2, 1.0)]

# ... which is not allowed
In [32]: try:
   ....:     ddemog = fwdpy11.DiscreteDemography(mass_migrations=move)
   ....: except ValueError as e:
   ....:     print(e)
   ....: 
DiscreteDemography: at time 0, attempting to move 200% of deme 0 is invalid

The rate of drift

Moving versus copying individuals is an important modeling choice. When you move individuals from one deme to another, the rate of drift changes in the source deme (as its size is reduced). This reduction in size is also a sudden bottleneck.

Copying, on the other hand, does not change the rate of drift in the source deme. However, it does seem to imply some sudden increase in fecundity that both came from nowhere and was short-lived.

Instantaneous deme size changes

Instantaneous changes in deme size are managed by instances of fwdpy11.SetDemeSize.

This class is relatively straightforward to use, so let’s dive right in:

In [33]: pop = fwdpy11.DiploidPopulation([20, 20], 1.0)

In [34]: dd = fwdpy11.DiscreteDemography(
   ....:     set_deme_sizes=[fwdpy11.SetDemeSize(when=5, deme=1, new_size=100)]
   ....: )
   ....: 

In [35]: st = SizeTracker()

In [36]: setup_and_run_model(pop, dd, 10, st)

In [37]: for i in st.data:
   ....:     print(i)
   ....: 
(1, 40, (array([0, 1], dtype=int32), array([20, 20])))
(2, 40, (array([0, 1], dtype=int32), array([20, 20])))
(3, 40, (array([0, 1], dtype=int32), array([20, 20])))
(4, 40, (array([0, 1], dtype=int32), array([20, 20])))
(5, 40, (array([0, 1], dtype=int32), array([20, 20])))
(6, 120, (array([0, 1], dtype=int32), array([ 20, 100])))
(7, 120, (array([0, 1], dtype=int32), array([ 20, 100])))
(8, 120, (array([0, 1], dtype=int32), array([ 20, 100])))
(9, 120, (array([0, 1], dtype=int32), array([ 20, 100])))
(10, 120, (array([0, 1], dtype=int32), array([ 20, 100])))

You may also kill off demes by setting their size to zero:

In [38]: pop = fwdpy11.DiploidPopulation([20, 20, 20], 1.0)

In [39]: dd = fwdpy11.DiscreteDemography(
   ....:     set_deme_sizes=[fwdpy11.SetDemeSize(when=5, deme=1, new_size=0)]
   ....: )
   ....: 

In [40]: st = SizeTracker()

In [41]: setup_and_run_model(pop, dd, 6, st)

In [42]: for i in st.data:
   ....:     print(i)
   ....: 
(1, 60, (array([0, 1, 2], dtype=int32), array([20, 20, 20])))
(2, 60, (array([0, 1, 2], dtype=int32), array([20, 20, 20])))
(3, 60, (array([0, 1, 2], dtype=int32), array([20, 20, 20])))
(4, 60, (array([0, 1, 2], dtype=int32), array([20, 20, 20])))
(5, 60, (array([0, 1, 2], dtype=int32), array([20, 20, 20])))
(6, 40, (array([0, 2], dtype=int32), array([20, 20])))

Changing growth rates

Instances of fwdpy11.SetExponentialGrowth manage the exponential growth rates per deme. Growth rates less than one indicate population decline, greater than one means growth and fwdpy11.NOGROWTH is equal to 1.0 to indicate no growth.

Let’s look at an example:

In [43]: pop = fwdpy11.DiploidPopulation([50], 1.0)

In [44]: g = [fwdpy11.SetExponentialGrowth(when=0, deme=0, G=1.1)]

In [45]: dd = fwdpy11.DiscreteDemography(set_growth_rates=g)

In [46]: st = SizeTracker()

In [47]: setup_and_run_model(pop, dd, 6, st)

In [48]: for i in st.data:
   ....:     print(i)
   ....: 
(1, 55, (array([0], dtype=int32), array([55])))
(2, 61, (array([0], dtype=int32), array([61])))
(3, 67, (array([0], dtype=int32), array([67])))
(4, 73, (array([0], dtype=int32), array([73])))
(5, 81, (array([0], dtype=int32), array([81])))
(6, 89, (array([0], dtype=int32), array([89])))

The deme sizes each generation must be integer values. The simulation uses C/C++ rules for rounding double-precision values to integer values. The function numpy.rint uses the same rules:

In [49]: N0 = np.float(50.0)

In [50]: for i in range(6):
   ....:     Ni = N0 * np.power(1.1, i + 1)
   ....:     print(i + 1, Ni, np.rint(Ni))
   ....: 
1 55.00000000000001 55.0
2 60.50000000000001 61.0
3 66.55000000000003 67.0
4 73.20500000000003 73.0
5 80.52550000000002 81.0
6 88.57805000000005 89.0

You may need to keep the rounding policy in mind when trying to predict final deme sizes when testing or when trying to convert a model from continuous time into discrete time.

Changing the selfing rate

Instances of fwdpy11.SetSelfingRate affect the rate of selfing-versus-outcrossing in different demes, or to change the rate within a deme over time. The default is that individuals don’t self unless they are picked twice as a parent by chance.

Using this type is straightforward. Before we dive in, we will create a new recorder type to track parents each generation:

# fmt: off
In [51]: class ParentTracker(object):
   ....:     def __init__(self):
   ....:         self.data = []
   ....:     def __call__(self, pop, sampler):
   ....:         for i in pop.diploid_metadata:
   ....:             self.data.append((i.label, i.deme, i.parents))
   ....: 

# fmt: on

Let’s run a simulation for a couple of generations:

In [52]: pop = fwdpy11.DiploidPopulation([5, 5], 1.0)

In [53]: sr = [fwdpy11.SetSelfingRate(when=0, deme=1, S=1.0)]  # Deme 1 always selfs

In [54]: dd = fwdpy11.DiscreteDemography(set_selfing_rates=sr)

In [55]: pt = ParentTracker()

In [56]: setup_and_run_model(pop, dd, 2, pt)

In our output, the deme label is the second value in each tuple, and any individual in deme 1 has the same parent listed twice because they were the product of a selfing event:

In [57]: for i in pt.data:
   ....:     print(i)
   ....: 
(0, 0, (1, 1))
(1, 0, (4, 1))
(2, 0, (0, 1))
(3, 0, (0, 4))
(4, 0, (1, 2))
(5, 1, (8, 8))
(6, 1, (9, 9))
(7, 1, (6, 6))
(8, 1, (5, 5))
(9, 1, (5, 5))
(0, 0, (0, 3))
(1, 0, (0, 2))
(2, 0, (2, 4))
(3, 0, (3, 3))
(4, 0, (0, 0))
(5, 1, (7, 7))
(6, 1, (8, 8))
(7, 1, (8, 8))
(8, 1, (5, 5))
(9, 1, (8, 8))

(In the above output, the parent IDs are the indexes of the parental individuals from their generation.)

Migration

For models with multiple demes, migration between then is managed by an instance of fwdpy11.MigrationMatrix.

For a migration matrix M, the default interpretation of M[i, j] is the fraction of deme i that will be replaced by migrations from deme j. The entry M[i, i] represents the non-migrant fraction of deme i’s ancestry. The matrix is “row-major” meaning that rows refer to migration into source demes. This definition of the migration matrix corresponds to that found in several different sources ([Christiansen1974], [Christiansen1975]). This definition of migration is also what diffusion models assume (e.g. [Jouganous2017]) as well as coalescent simulations like msprime [Kelleher2016].

For example, consider the following matrix:

In [58]: m = np.array([0.9, 0.1, 0.5, 0.5]).reshape(2, 2)

In [59]: m
Out[59]: 
array([[0.9, 0.1],
       [0.5, 0.5]])

The first row corresponds to the ancestry of deme 0, such that 90% of offspring will be non-migrants and 10% will be migrants from deme 1:

In [60]: m[
   ....:     0,
   ....: ]
   ....: 
Out[60]: array([0.9, 0.1])

To be concrete, if the size of deme 0 in the next generation is 1,000, then the expected number of migrant and non-migrant offspring of offspring in deme 0 is:

In [61]: m[0,] * 1e3
Out[61]: array([900., 100.])

The second row implies that half the ancestry of deme 1 is due to migrants and half due to non-migrants:

In [62]: m[
   ....:     1,
   ....: ]
   ....: 
Out[62]: array([0.5, 0.5])

The numpy array is sufficient to construct our demographic model:

In [63]: d = fwdpy11.DiscreteDemography(migmatrix=m)

In [64]: d.migmatrix
Out[64]: 
fwdpy11.MigrationMatrix(migmatrix=array([[0.9, 0.1],
       [0.5, 0.5]]), scaled=False)

By default, there is no migration, which is represented by the value None. For example, the following model has no migration events:

# Define demographic events w/o any migration stuff
In [65]: d = fwdpy11.DiscreteDemography(set_deme_sizes=[fwdpy11.SetDemeSize(0, 1, 500)])

In [66]: d.migmatrix is None
Out[66]: True

In order to specify a model with no initial migration, you may use an identity matrix:

In [67]: d = fwdpy11.DiscreteDemography(migmatrix=np.identity(2))

In [68]: d.migmatrix
Out[68]: 
fwdpy11.MigrationMatrix(migmatrix=array([[1., 0.],
       [0., 1.]]), scaled=False)

The only reason to use the identity matrix is to start a simulation with no migration and then change the rates later via instances of fwdpy11.SetMigrationRates. To see this in action, we’ll first generate a new type to track if parents of offspring in deme 1 are migrants or not:

# fmt: off
In [69]: class MigrationTracker(object):
   ....:     def __init__(self, N0):
   ....:         self.N0 = N0
   ....:         self.data = []
   ....:     def __call__(self, pop, sampler):
   ....:         for i in pop.diploid_metadata:
   ....:             if i.deme == 1:
   ....:                 p = []
   ....:                 for j in i.parents:
   ....:                     if j < self.N0:
   ....:                         p.append((j, True))
   ....:                     else:
   ....:                         p.append((j, False))
   ....:                 self.data.append((pop.generation, i.label, i.deme, p))
   ....: 

# fmt: on
# No migration at first
In [70]: mm = np.identity(2)

# In generation 3, reset migration rates for deme 1 such
# that parents are equally likey from both demes.
In [71]: cm = [fwdpy11.SetMigrationRates(3, 1, [0.5, 0.5])]

In [72]: dd = fwdpy11.DiscreteDemography(migmatrix=mm, set_migration_rates=cm)

In [73]: pop = fwdpy11.DiploidPopulation([10, 10], 1.0)

In [74]: mt = MigrationTracker(10)

In [75]: setup_and_run_model(pop, dd, 4, mt)
In [76]: for i in mt.data:
   ....:     nmig = 0
   ....:     if i[1] > 10:
   ....:         if i[3][0][1] is True:
   ....:             nmig += 1
   ....:         if i[3][1][1] is True:
   ....:             nmig += 1
   ....:     mstring = ""
   ....:     if nmig > 0:
   ....:         mstring = "<- {} migrant parent".format(nmig)
   ....:     if nmig > 1:
   ....:         mstring += "s"
   ....:     print(i, mstring)
   ....: 
(1, 10, 1, [(17, False), (15, False)]) 
(1, 11, 1, [(14, False), (14, False)]) 
(1, 12, 1, [(15, False), (15, False)]) 
(1, 13, 1, [(10, False), (19, False)]) 
(1, 14, 1, [(17, False), (17, False)]) 
(1, 15, 1, [(16, False), (10, False)]) 
(1, 16, 1, [(16, False), (10, False)]) 
(1, 17, 1, [(19, False), (14, False)]) 
(1, 18, 1, [(17, False), (10, False)]) 
(1, 19, 1, [(17, False), (10, False)]) 
(2, 10, 1, [(16, False), (13, False)]) 
(2, 11, 1, [(12, False), (17, False)]) 
(2, 12, 1, [(12, False), (10, False)]) 
(2, 13, 1, [(13, False), (11, False)]) 
(2, 14, 1, [(15, False), (14, False)]) 
(2, 15, 1, [(10, False), (19, False)]) 
(2, 16, 1, [(15, False), (10, False)]) 
(2, 17, 1, [(15, False), (15, False)]) 
(2, 18, 1, [(12, False), (12, False)]) 
(2, 19, 1, [(15, False), (15, False)]) 
(3, 10, 1, [(15, False), (16, False)]) 
(3, 11, 1, [(18, False), (17, False)]) 
(3, 12, 1, [(15, False), (11, False)]) 
(3, 13, 1, [(12, False), (10, False)]) 
(3, 14, 1, [(12, False), (18, False)]) 
(3, 15, 1, [(14, False), (10, False)]) 
(3, 16, 1, [(15, False), (11, False)]) 
(3, 17, 1, [(12, False), (16, False)]) 
(3, 18, 1, [(13, False), (10, False)]) 
(3, 19, 1, [(19, False), (16, False)]) 
(4, 10, 1, [(0, True), (5, True)]) 
(4, 11, 1, [(1, True), (4, True)]) <- 2 migrant parents
(4, 12, 1, [(0, True), (7, True)]) <- 2 migrant parents
(4, 13, 1, [(5, True), (2, True)]) <- 2 migrant parents
(4, 14, 1, [(0, True), (6, True)]) <- 2 migrant parents
(4, 15, 1, [(14, False), (14, False)]) 
(4, 16, 1, [(18, False), (12, False)]) 
(4, 17, 1, [(11, False), (11, False)]) 
(4, 18, 1, [(0, True), (8, True)]) <- 2 migrant parents
(4, 19, 1, [(11, False), (18, False)]) 

An alternative model of migration

The description of migration rates above implies that migration events are independent of of source deme sizes. To revisit our earlier example:

In [77]: m = np.array([0.9, 0.1, 0.5, 0.5]).reshape(2, 2)

# The is the expected number of parents from demes 0 and 1
# to offspring born in deme 0:
In [78]: m[0,] * 1000
Out[78]: array([900., 100.])

fwdpy11 allows for a different migration scheme where the size of the source deme matters. For this model, M[i ,j] is the probability that an individual with parents from deme j is born in deme i. Internally, the migration matrix entries M[i, j] are multiplied by the size of the source demes, which means that larger demes with nonzero migration rates to other demes have a larger chance of being sources of migrant offspring.

For example:

In [79]: deme_sizes = np.array([1000, 2000])

In [80]: m
Out[80]: 
array([[0.9, 0.1],
       [0.5, 0.5]])

In [81]: md = m * deme_sizes

# The following line divides each
# row by its sum
In [82]: md / np.sum(md, axis=1)[:, None]
Out[82]: 
array([[0.81818182, 0.18181818],
       [0.33333333, 0.66666667]])

The first matrix is the same as in the preceding section–90% of the offspring in deme 0 will have parents from deme 0. In the second matrix, that fraction is reduced to about 82% because deme 1 is twice as large as deme 0.

To enable this migration model, create an instance of fwdpy11.MigrationMatrix and pass True as the second parameter:

In [83]: M = fwdpy11.MigrationMatrix(m, True)

In [84]: d = fwdpy11.DiscreteDemography(migmatrix=M)

This will also work, but is less explicit:

In [85]: d = fwdpy11.DiscreteDemography(migmatrix=(m, True))

Note

This model of migration will typically give different results from diffusion models and coalescent simulations!

Examples of models

Isolation with migration, or “IM”

Consider two demes that split apart T time units ago and then grow to different sizes in the present. After the split, migration occurs between the two demes. The demographic model has the following parameters:

  • Nanc, the ancestral population size.

  • T, the time of the split, which is in units of Nanc.

  • psplit, the proportion of the ancestral population that splits off to found deme 1.

  • N0, the final size of deme 0, relative to Nanc.

  • N1, the final size of deme 1, relative to Nanc.

  • m01, the migration rate from deme 0 to deme 1.

  • m10, the migration rate from deme 1 to deme 0.

Here is the model in its entirety, with no mutation and no recombination. First, we will set up the demographic events. The population with evolve for Nanc generations before the split.

In [86]: Nanc = 100

In [87]: T = 0.2

In [88]: psplit = 0.33

In [89]: N0, N1 = 2, 3

In [90]: m01, m10 = 0.01, 0.0267

# The split event
In [91]: split = [fwdpy11.move_individuals(when=Nanc, source=0, destination=1, fraction=psplit)]

# Get growth rates and set growth rate changes,
# taking care to handle our rounding!
In [92]: gens_post_split = np.rint(Nanc * T).astype(int)

In [93]: N0split = np.rint(Nanc * (1.0 - psplit))

In [94]: N0final = np.rint(N0 * Nanc)

In [95]: N1split = np.rint(Nanc * psplit)

In [96]: N1final = np.rint(N1 * Nanc)

In [97]: G0 = fwdpy11.exponential_growth_rate(N0split, N0final, gens_post_split)

In [98]: G1 = fwdpy11.exponential_growth_rate(N1split, N1final, gens_post_split)

In [99]: growth = [
   ....:     fwdpy11.SetExponentialGrowth(Nanc, 0, G0),
   ....:     fwdpy11.SetExponentialGrowth(Nanc, 1, G1),
   ....: ]
   ....: 

# Set up the migration matrix for two demes, but only
# deme zero exists.
In [100]: m = fwdpy11.migration_matrix_single_extant_deme(2, 0)

# The rows of the matrix change at the split:
In [101]: cm = [
   .....:     fwdpy11.SetMigrationRates(Nanc, 0, [1.0 - m10, m10]),
   .....:     fwdpy11.SetMigrationRates(Nanc, 1, [m01, 1.0 - m01]),
   .....: ]
   .....: 

In [102]: d = fwdpy11.DiscreteDemography(
   .....:     mass_migrations=split, set_growth_rates=growth, set_migration_rates=cm, migmatrix=m
   .....: )
   .....: 

The above code made use of two helper functions:

Finally, we can run it:

In [103]: pop = fwdpy11.DiploidPopulation(Nanc, 1.0)

In [104]: setup_and_run_model(pop, d, Nanc + gens_post_split)

Now we check the final population sizes and make sure they are correct:

In [105]: ds = pop.deme_sizes()

In [106]: assert ds[1][0] == N0final

In [107]: assert ds[1][1] == N1final

This model is common enough that you shouldn’t have to implement it from scratch each time. For this reason, we provide it in fwdpy11.demographic_models.IM.two_deme_IM().

In [108]: import fwdpy11.demographic_models.IM

In [109]: dmodel = fwdpy11.demographic_models.IM.two_deme_IM(
   .....:     Nanc, T, psplit, (N0, N1), (m01, m10), burnin=1.0
   .....: )
   .....: 

In [110]: pop2 = fwdpy11.DiploidPopulation(Nanc, 1.0)

In [111]: setup_and_run_model(pop2, dmodel.model, dmodel.metadata.simlen)

In [112]: assert pop.generation == pop2.generation

In [113]: assert pop2.generation == dmodel.metadata.simlen

In [114]: ds2 = pop2.deme_sizes()

In [115]: assert np.array_equal(ds[0], ds2[0])

In [116]: assert np.array_equal(ds[1], ds2[1])

See Isolation-with-migration (IM) for an example of using this function to compare to results from diffusion models.