Skip to content

Commit

Permalink
reimplemented AncientEurope model with msprime demography API
Browse files Browse the repository at this point in the history
  • Loading branch information
@awohns @mufernando
awohns authored and mufernando committed May 19, 2022
1 parent bfded34 commit e6885e8
Showing 1 changed file with 78 additions and 74 deletions.
152 changes: 78 additions & 74 deletions stdpopsim/catalog/HomSap/demographic_models.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1820,15 +1820,22 @@ def _ancient_europe():
from each population.
"""
populations = [
stdpopsim.Population(id="OOA", description="Basal/OOA"),
stdpopsim.Population(id="NE", description="Northern European"),
stdpopsim.Population(id="WA", description="West Asian"),
stdpopsim.Population(
id="Pop0",
description="1000GenomesEUR/BronzeAge/Neolithic/Anatolian/WestAsian/Basal",
id="CHG", description="Caucasus Hunter-gathers", sampling_time=300
),
stdpopsim.Population(id="Pop1", description="Yamnaya/CHG", sampling_time=140),
stdpopsim.Population(id="ANA", description="Anatolian", sampling_time=260),
stdpopsim.Population(
id="Pop2", description="WHG/NorthernEuropean", sampling_time=200
id="WHG", description="Western Hunter-gathers", sampling_time=250
),
stdpopsim.Population(id="Pop3", description="EHG", sampling_time=180),
stdpopsim.Population(
id="EHG", description="Eastern Hunter-gathers", sampling_time=250
),
stdpopsim.Population(id="NEO", description="Neolithic", sampling_time=180),
stdpopsim.Population(id="YAM", description="Yamnaya", sampling_time=160),
stdpopsim.Population(id="Bronze", description="Bronze Age", sampling_time=135),
]
citations = [
stdpopsim.Citation(
Expand All @@ -1843,95 +1850,92 @@ def _ancient_europe():
mutation_rate = 1.25e-8

# initial population sizes:
N_bronze = 50000
N_Yam = 5000
N_whg = 10000
N_ehg = 10000
N_neo = 50000
N_chg = 10000
N_Bronze = 50000
N_YAM = 5000
N_WHG = 10000
N_EHG = 10000
N_ANA = 50000
N_NEO = 50000
N_CHG = 10000
N_NE = 5000 # Ancestor of WHG and EHG
N_WA = 5000 # Ancestor of CHG and Anatolian farmers
N_OOA = 5000

# Time of events
T_bronze = 140
T_Yam = 180
T_neo = 200
T_near_east = 800
T_europe = 600
T_basal = 1500
T_Bronze = 140
T_YAM = 180
T_NEO = 200
T_Near_East = 800
T_Europe = 600
T_Basal = 1500

# Growth rate and initial population size for present day from bronze age
r_EU = 0.067
N_present = N_bronze / math.exp(-r_EU * T_bronze)
N_Present = N_Bronze / math.exp(-r_EU * T_Bronze)

population_configurations = [
msprime.PopulationConfiguration(
initial_size=N_present, growth_rate=r_EU, metadata=populations[0].asdict()
),
msprime.PopulationConfiguration(
initial_size=N_Yam, metadata=populations[1].asdict()
),
msprime.PopulationConfiguration(
initial_size=N_whg, metadata=populations[2].asdict()
),
msprime.PopulationConfiguration(
initial_size=N_ehg, metadata=populations[3].asdict()
),
]
# Create msprime demography object
de = msprime.Demography()

Bronze_formation = [
msprime.MassMigration(time=T_bronze, source=0, dest=1, proportion=0.5),
msprime.PopulationParametersChange(
time=T_bronze, initial_size=N_neo, growth_rate=0, population=0
),
]

Yam_formation = [
msprime.MassMigration(time=T_Yam, source=1, dest=3, proportion=0.5),
msprime.PopulationParametersChange(
time=T_Yam, initial_size=N_chg, population=1
),
]

European_neolithic = [
msprime.MassMigration(time=T_neo, source=0, dest=2, proportion=1.0 / 4.0)
]

HG_split = [
msprime.MassMigration(time=T_europe, source=3, dest=2, proportion=1),
msprime.PopulationParametersChange(
time=T_europe, initial_size=N_NE, population=2
),
]

Near_east_split = [
msprime.MassMigration(time=T_near_east, source=1, dest=0, proportion=1),
msprime.PopulationParametersChange(
time=T_near_east, initial_size=N_WA, population=0
),
]
# Create populations
de.add_population(
initial_size=N_OOA, name="OOA", description=populations[0].description
)
de.add_population(
initial_size=N_NE, name="NE", description=populations[1].description
)
de.add_population(
initial_size=N_WA, name="WA", description=populations[2].description
)
de.add_population(
initial_size=N_CHG, name="CHG", description=populations[3].description
)
de.add_population(
initial_size=N_ANA, name="ANA", description=populations[4].description
)
de.add_population(
initial_size=N_WHG, name="WHG", description=populations[5].description
)
de.add_population(
initial_size=N_EHG, name="EHG", description=populations[6].description
)
de.add_population(
initial_size=N_YAM, name="YAM", description=populations[7].description
)

Basal_split = [msprime.MassMigration(time=T_basal, source=2, dest=0, proportion=1)]
de.add_population(
initial_size=N_NEO, name="NEO", description=populations[8].description
)
de.add_population(
initial_size=N_Present,
growth_rate=r_EU,
name="Bronze",
description=populations[9].description,
)

demographic_events = (
Bronze_formation
+ Yam_formation
+ European_neolithic
+ HG_split
+ Near_east_split
+ Basal_split
de.add_admixture(
time=T_Bronze,
derived="Bronze",
ancestral=["YAM", "NEO"],
proportions=[0.5, 0.5],
)
de.add_admixture(
time=T_YAM, derived="YAM", ancestral=["EHG", "CHG"], proportions=[0.5, 0.5]
)
de.add_admixture(
time=T_NEO, derived="NEO", ancestral=["WHG", "ANA"], proportions=[0.25, 0.75]
)
de.add_population_split(time=T_Europe, derived=["WHG", "EHG"], ancestral="NE")
de.add_population_split(time=T_Near_East, derived=["CHG", "ANA"], ancestral="WA")
de.add_population_split(time=T_Basal, derived=["NE", "WA"], ancestral="OOA")

return stdpopsim.DemographicModel(
id=id,
description=description,
long_description=long_description,
populations=populations,
model=de,
citations=citations,
generation_time=generation_time,
mutation_rate=mutation_rate,
population_configurations=population_configurations,
demographic_events=demographic_events,
)


Expand Down

0 comments on commit e6885e8

Please sign in to comment.