model.py

import math

#start with population of size 10k
N_AF = 10000
#when we split this into five, we should have initial size before growth of N/5 in each
N_EU0 = N_AF / 5

#the five populations merged into the ancestral population 100 generations ago
T_merge = 100
#assume growth rate of 0.4%
r_EU = 0.004
#calculate present day size based on starting size, growth time and growth rate 
N_EU = N_EU0 / math.exp(-r_EU * T_merge)

#migration rate between populations
mEU = 9.6e-5

population_configurations = [
    #0: ancestral population - we don't sample from this
    msprime.PopulationConfiguration(sample_size=0, initial_size=N_AF, growth_rate=0),
    
    #1-5: modeled after EUR (initial size = current size, after growth, bottom of tree)
    msprime.PopulationConfiguration(sample_size=5000, initial_size=N_EU, growth_rate=r_EU),
    msprime.PopulationConfiguration(sample_size=5000, initial_size=N_EU, growth_rate=r_EU),
    msprime.PopulationConfiguration(sample_size=5000, initial_size=N_EU, growth_rate=r_EU),
    msprime.PopulationConfiguration(sample_size=5000, initial_size=N_EU, growth_rate=r_EU),
    msprime.PopulationConfiguration(sample_size=5000, initial_size=N_EU, growth_rate=r_EU),
]
migration_matrix = [
    #no migration from/to ancestral population, mEU migration between the others
    [0, 0, 0, 0, 0, 0],
    [0, 0, mEU, mEU, mEU, mEU],
    [0, mEU, 0, mEU, mEU, mEU],
    [0, mEU, mEU, 0, mEU, mEU],
    [0, mEU, mEU, mEU, 0, mEU],
    [0, mEU, mEU, mEU, mEU, 0],
]
demographic_events = [
    #merge the populations into the ancestral population
    msprime.MassMigration(time=T_merge, source=1, destination=0, proportion=1.0),
    msprime.MassMigration(time=T_merge, source=2, destination=0, proportion=1.0),
    msprime.MassMigration(time=T_merge, source=3, destination=0, proportion=1.0),
    msprime.MassMigration(time=T_merge, source=4, destination=0, proportion=1.0),
    msprime.MassMigration(time=T_merge, source=5, destination=0, proportion=1.0),
]
# Use the demography debugger to print out the demographic history
# that we have just described.
dp = msprime.DemographyDebugger(
    Ne=N_AF,
    population_configurations=population_configurations,
    migration_matrix=migration_matrix,
    demographic_events=demographic_events)
dp.print_history()

"""
=============================
Epoch: 0 -- 100.0 generations
=============================
     start     end      growth_rate |     0        1        2        3        4        5    
   -------- --------       -------- | -------- -------- -------- -------- -------- -------- 
0 |  1e+04    1e+04               0 |     0        0        0        0        0        0    
1 |2.98e+03   2e+03           0.004 |     0        0     9.6e-05  9.6e-05  9.6e-05  9.6e-05 
2 |2.98e+03   2e+03           0.004 |     0     9.6e-05     0     9.6e-05  9.6e-05  9.6e-05 
3 |2.98e+03   2e+03           0.004 |     0     9.6e-05  9.6e-05     0     9.6e-05  9.6e-05 
4 |2.98e+03   2e+03           0.004 |     0     9.6e-05  9.6e-05  9.6e-05     0     9.6e-05 
5 |2.98e+03   2e+03           0.004 |     0     9.6e-05  9.6e-05  9.6e-05  9.6e-05     0    

Events @ generation 100.0
   - Mass migration: lineages move from 1 to 0 with probability 1.0
   - Mass migration: lineages move from 2 to 0 with probability 1.0
   - Mass migration: lineages move from 3 to 0 with probability 1.0
   - Mass migration: lineages move from 4 to 0 with probability 1.0
   - Mass migration: lineages move from 5 to 0 with probability 1.0

===============================
Epoch: 100.0 -- inf generations
===============================
     start     end      growth_rate |     0        1        2        3        4        5    
   -------- --------       -------- | -------- -------- -------- -------- -------- -------- 
0 |  1e+04    1e+04               0 |     0        0        0        0        0        0    
1 |2.98e+03     0             0.004 |     0        0     9.6e-05  9.6e-05  9.6e-05  9.6e-05 
2 |2.98e+03     0             0.004 |     0     9.6e-05     0     9.6e-05  9.6e-05  9.6e-05 
3 |2.98e+03     0             0.004 |     0     9.6e-05  9.6e-05     0     9.6e-05  9.6e-05 
4 |2.98e+03     0             0.004 |     0     9.6e-05  9.6e-05  9.6e-05     0     9.6e-05 
5 |2.98e+03     0             0.004 |     0     9.6e-05  9.6e-05  9.6e-05  9.6e-05     0    
"""