# 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:

Start with 50 individuals and

*copy*them to deme 1 in generation 0Start 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.