Repository: popsim-consortium/stdpopsim (language: python)

Issue:
Refactor catalog code into submodules
The ``__init__.py`` file in the catalog subpackages currently contains a lot of code, which will confuse contributors I think. It would be better if this was split into separate files, .e.g., ``__init__.py`` would become

```
from . import genome
from . import genetic_maps
from . import annotations
from . import demographic_models
```
This would probably help contributors see where their code is supposed to go and make things less confusing.

Thoughts @grahamgower? It's not hard to do, just a bit of file rejiggering.


Discussion / hints:
Please do! That's much nicer!

Fix patch:
diff --git a/stdpopsim/catalog/AraTha/__init__.py b/stdpopsim/catalog/AraTha/__init__.py
index 883cfa19f..977e5cc0c 100644
--- a/stdpopsim/catalog/AraTha/__init__.py
+++ b/stdpopsim/catalog/AraTha/__init__.py
@@ -1,366 +1,6 @@
 """
-Genome, genetic map and demographic model definitions for Arabidopsis thaliana.
+Catalog definitions for Arabidopsis Thaliana.
 """
-import msprime
-import numpy as np
-
-import stdpopsim
-from . import genome_data
-
-###########################################################
-#
-# Genome definition
-#
-###########################################################
-
-_mean_recombination_rate = 200 / 124000 / 2 / 1e6
-
-_recombination_rate_data = {str(j): _mean_recombination_rate for j in range(1, 6)}
-_recombination_rate_data["Mt"] = 0
-_recombination_rate_data["Pt"] = 0  # JK Is this correct??
-
-
-# mutation rate from Ossowski 2010 Science
-# recombination value from Huber et al 2014 MBE
-# rho=200/Mb, assume Ne=124,000, rho=2*Ne*r
-_chromosomes = []
-for name, data in genome_data.data["chromosomes"].items():
-    _chromosomes.append(
-        stdpopsim.Chromosome(
-            id=name,
-            length=data["length"],
-            synonyms=data["synonyms"],
-            mutation_rate=7e-9,
-            recombination_rate=_recombination_rate_data[name],
-        )
-    )
-
-_genome = stdpopsim.Genome(
-    chromosomes=_chromosomes,
-    assembly_name=genome_data.data["assembly_name"],
-    assembly_accession=genome_data.data["assembly_accession"],
-    mutation_rate_citations=[
-        stdpopsim.Citation(
-            author="Ossowski et al.",
-            year="2010",
-            doi="https://doi.org/10.1126/science.1180677",
-            reasons={stdpopsim.CiteReason.MUT_RATE},
-        )
-    ],
-    recombination_rate_citations=[
-        stdpopsim.Citation(
-            author="Huber et al.",
-            year="2014",
-            doi="https://doi.org/10.1093/molbev/msu247",
-            reasons={stdpopsim.CiteReason.REC_RATE},
-        )
-    ],
-    assembly_citations=[
-        stdpopsim.Citation(
-            doi="https://doi.org/10.1093/nar/gkm965",
-            year="2007",
-            author="Swarbreck et al.",
-            reasons={stdpopsim.CiteReason.ASSEMBLY},
-        )
-    ],
-)
-
-_species = stdpopsim.Species(
-    id="AraTha",
-    name="Arabidopsis thaliana",
-    common_name="A. thaliana",
-    genome=_genome,
-    generation_time=1.0,
-    generation_time_citations=[
-        stdpopsim.Citation(
-            doi="https://doi.org/10.1890/0012-9658(2002)083[1006:GTINSO]2.0.CO;2",
-            year="2002",
-            author="Donohue",
-            reasons={stdpopsim.CiteReason.GEN_TIME},
-        )
-    ],
-    population_size=10 ** 4,
-    population_size_citations=[
-        stdpopsim.Citation(
-            doi="https://doi.org/10.1016/j.cell.2016.05.063",
-            year="2016",
-            author="1001GenomesConsortium",
-            reasons={stdpopsim.CiteReason.POP_SIZE},
-        )
-    ],
-)
-
-stdpopsim.register_species(_species)
-
-###########################################################
-#
-# Genetic maps
-#
-###########################################################
-
-_gm = stdpopsim.GeneticMap(
-    species=_species,
-    id="SalomeAveraged_TAIR7",
-    description="Crossover frequency map averaged over 17 populations",
-    long_description="""
-        This map is based on the study of crossover frequencies in over 7000
-        plants in 17 F2 populations derived from crosses between 18 A. thaliana
-        accessions. Salomé et al provide genetic maps for each of these
-        populations. To get a single map for each chromosome, the Haldane map
-        function distances were converted to recombination rates (cM/Mb) for
-        each cross and then averaged across the 17 populations using loess.
-        """,
-    url=(
-        "https://stdpopsim.s3-us-west-2.amazonaws.com/genetic_maps/"
-        "AraTha/salome2012_maps.tar.gz"
-    ),
-    sha256="49745e1cab87d59e33eacfdf66303839632d3b07883dd55a99fe1dc27b336ac6",
-    file_pattern="arab_chr{id}_map_loess.txt",
-    citations=[
-        stdpopsim.Citation(
-            doi="https://doi.org/10.1038/hdy.2011.95",
-            author="Salomé et al.",
-            year=2012,
-            reasons={stdpopsim.CiteReason.GEN_MAP},
-        )
-    ],
-)
-_species.add_genetic_map(_gm)
-
-
-###########################################################
-#
-# Demographic models
-#
-###########################################################
-
-
-def _sma_1pop():
-    # the size during the interval times[k] to times[k+1] = sizes[k]
-    times = np.array(
-        [
-            699,
-            2796,
-            6068,
-            9894,
-            14370,
-            19606,
-            25730,
-            32894,
-            41275,
-            51077,
-            62544,
-            75958,
-            91648,
-            110001,
-            131471,
-            156584,
-            185960,
-            220324,
-            260520,
-            307540,
-            362541,
-            426879,
-            502139,
-            590173,
-            693151,
-            813610,
-            954517,
-            1119341,
-            1312147,
-            1537686,
-            1801500,
-            2110100,
-        ]
-    )
-    sizes = np.array(
-        [
-            42252426,
-            42252426,
-            60323,
-            72174,
-            40591,
-            21158,
-            21442,
-            39942,
-            78908,
-            111132,
-            110745,
-            96283,
-            87661,
-            83932,
-            83829,
-            91813,
-            111644,
-            143456,
-            181571,
-            217331,
-            241400,
-            246984,
-            238593,
-            228222,
-            217752,
-            198019,
-            165210,
-            121796,
-            121796,
-            73989,
-            73989,
-            73989,
-        ]
-    )
-
-    # MSMC is accurate from 40Kya-1.6Mya for A.thaliana (Durvasula et al 2017)
-    # set the first 7 sizes
-    # equal to the size at 8 (~40Kya)
-    sizes[:8] = sizes[8]
-    # set the last 2 entries equal
-    # to the size at 30 (~1.6Mya)
-    sizes[30:32] = sizes[30]
-
-    demographic_events = []
-    for sz, t in zip(sizes, times):
-        demographic_events.append(
-            msprime.PopulationParametersChange(time=t, initial_size=sz, population_id=0)
-        )
-
-    populations = [
-        stdpopsim.Population(
-            id="SouthMiddleAtlas",
-            description="Arabidopsis Thaliana South Middle Atlas population",
-        )
-    ]
-
-    return stdpopsim.DemographicModel(
-        id="SouthMiddleAtlas_1D17",
-        description="South Middle Atlas piecewise constant size",
-        long_description="""
-            This model comes from MSMC using two randomly sampled homozygous
-            individuals (Khe32 and Ifr4) from the South Middle Atlas region
-            from the Middle Atlas Mountains in Morocco. The model is estimated
-            with 32 time periods. Because estimates from the recent and ancient
-            past are less accurate, we set the population size in the first 7
-            time periods equal to the size at the 8th time period and the size
-            during last 2 time periods equal to the size in the 30th time
-            period.
-        """,
-        populations=populations,
-        citations=[
-            stdpopsim.Citation(
-                author="Durvasula et al.",
-                year=2017,
-                doi="https://doi.org/10.1073/pnas.1616736114",
-                reasons={stdpopsim.CiteReason.DEM_MODEL},
-            )
-        ],
-        generation_time=1,
-        demographic_events=demographic_events,
-        population_configurations=[
-            msprime.PopulationConfiguration(
-                initial_size=sizes[0], metadata=populations[0].asdict()
-            )
-        ],
-    )
-
-
-_species.add_demographic_model(_sma_1pop())
-
-
-def _afr_2epoch():
-    N_A = 746148
-    N_0 = 100218
-    t_1 = 568344
-    populations = [
-        stdpopsim.Population(
-            id="SouthMiddleAtlas",
-            description="Arabidopsis Thaliana South Middle Atlas population",
-        )
-    ]
-    return stdpopsim.DemographicModel(
-        id="African2Epoch_1H18",
-        description="South Middle Atlas African two epoch model",
-        long_description="""
-            Model estimated from site frequency spectrum of synonymous
-            SNPs from African South Middle Atlas samples using
-            Williamson et al. 2005 methodology. Values come from supplementary
-            table 1 of Huber et al 2018. Sizes change from N_A -> N_0 and t_1 is
-            time of the second epoch.
-        """,
-        populations=populations,
-        citations=[
-            stdpopsim.Citation(
-                author="Huber et al.",
-                year=2018,
-                doi="https://doi.org/10.1038/s41467-018-05281-7",
-                reasons={stdpopsim.CiteReason.DEM_MODEL},
-            )
-        ],
-        generation_time=1,
-        population_configurations=[
-            msprime.PopulationConfiguration(
-                initial_size=N_0, metadata=populations[0].asdict()
-            )
-        ],
-        demographic_events=[
-            msprime.PopulationParametersChange(
-                time=t_1, initial_size=N_A, population_id=0
-            )
-        ],
-    )
-
-
-_species.add_demographic_model(_afr_2epoch())
-
-
-def _afr_3epoch():
-    # values from the supplementary Table 1.
-    # the size changed as N_A -> N_2 -> N_3.
-    # t_2 is time of 2nd epoch and t_3 of the third epoch
-    N_A = 161744
-    N_2 = 24076
-    N_3 = 203077
-    t_2 = 7420
-    t_3 = 14534
-    populations = [
-        stdpopsim.Population(
-            id="SouthMiddleAtlas",
-            description="Arabidopsis Thaliana South Middle Atlas population",
-        )
-    ]
-    return stdpopsim.DemographicModel(
-        id="African3Epoch_1H18",
-        description="South Middle Atlas African three epoch model",
-        long_description="""
-            Model estimated from site frequency spectrum of synonymous
-            SNPs from African (South Middle Atlas) samples using Williamson et
-            al. 2005 methodology. Values come from supplementary table 1 of
-            Huber et al 2018. Sizes change from N_A -> N_2 -> N_3 and t_2 is
-            the time of the second epoch and t_3 is the time of the 3rd epoch.
-        """,
-        populations=populations,
-        citations=[
-            stdpopsim.Citation(
-                author="Huber et al.",
-                year=2018,
-                doi="https://doi.org/10.1038/s41467-018-05281-7",
-                reasons={stdpopsim.CiteReason.DEM_MODEL},
-            )
-        ],
-        generation_time=1,
-        population_configurations=[
-            msprime.PopulationConfiguration(
-                initial_size=N_3, metadata=populations[0].asdict()
-            )
-        ],
-        demographic_events=[
-            msprime.PopulationParametersChange(
-                time=t_3, initial_size=N_2, population_id=0
-            ),
-            msprime.PopulationParametersChange(
-                time=t_2 + t_3, initial_size=N_A, population_id=0
-            ),
-        ],
-    )
-
-
-_species.add_demographic_model(_afr_3epoch())
+from . import species  # noqa 401
+from . import genetic_maps  # noqa 401
+from . import demographic_models  # noqa 401
diff --git a/stdpopsim/catalog/AraTha/demographic_models.py b/stdpopsim/catalog/AraTha/demographic_models.py
new file mode 100644
index 000000000..037dc8173
--- /dev/null
+++ b/stdpopsim/catalog/AraTha/demographic_models.py
@@ -0,0 +1,237 @@
+import msprime
+import numpy as np
+
+import stdpopsim
+
+_species = stdpopsim.get_species("AraTha")
+
+
+def _sma_1pop():
+    # the size during the interval times[k] to times[k+1] = sizes[k]
+    times = np.array(
+        [
+            699,
+            2796,
+            6068,
+            9894,
+            14370,
+            19606,
+            25730,
+            32894,
+            41275,
+            51077,
+            62544,
+            75958,
+            91648,
+            110001,
+            131471,
+            156584,
+            185960,
+            220324,
+            260520,
+            307540,
+            362541,
+            426879,
+            502139,
+            590173,
+            693151,
+            813610,
+            954517,
+            1119341,
+            1312147,
+            1537686,
+            1801500,
+            2110100,
+        ]
+    )
+    sizes = np.array(
+        [
+            42252426,
+            42252426,
+            60323,
+            72174,
+            40591,
+            21158,
+            21442,
+            39942,
+            78908,
+            111132,
+            110745,
+            96283,
+            87661,
+            83932,
+            83829,
+            91813,
+            111644,
+            143456,
+            181571,
+            217331,
+            241400,
+            246984,
+            238593,
+            228222,
+            217752,
+            198019,
+            165210,
+            121796,
+            121796,
+            73989,
+            73989,
+            73989,
+        ]
+    )
+
+    # MSMC is accurate from 40Kya-1.6Mya for A.thaliana (Durvasula et al 2017)
+    # set the first 7 sizes
+    # equal to the size at 8 (~40Kya)
+    sizes[:8] = sizes[8]
+    # set the last 2 entries equal
+    # to the size at 30 (~1.6Mya)
+    sizes[30:32] = sizes[30]
+
+    demographic_events = []
+    for sz, t in zip(sizes, times):
+        demographic_events.append(
+            msprime.PopulationParametersChange(time=t, initial_size=sz, population_id=0)
+        )
+
+    populations = [
+        stdpopsim.Population(
+            id="SouthMiddleAtlas",
+            description="Arabidopsis Thaliana South Middle Atlas population",
+        )
+    ]
+
+    return stdpopsim.DemographicModel(
+        id="SouthMiddleAtlas_1D17",
+        description="South Middle Atlas piecewise constant size",
+        long_description="""
+            This model comes from MSMC using two randomly sampled homozygous
+            individuals (Khe32 and Ifr4) from the South Middle Atlas region
+            from the Middle Atlas Mountains in Morocco. The model is estimated
+            with 32 time periods. Because estimates from the recent and ancient
+            past are less accurate, we set the population size in the first 7
+            time periods equal to the size at the 8th time period and the size
+            during last 2 time periods equal to the size in the 30th time
+            period.
+        """,
+        populations=populations,
+        citations=[
+            stdpopsim.Citation(
+                author="Durvasula et al.",
+                year=2017,
+                doi="https://doi.org/10.1073/pnas.1616736114",
+                reasons={stdpopsim.CiteReason.DEM_MODEL},
+            )
+        ],
+        generation_time=1,
+        demographic_events=demographic_events,
+        population_configurations=[
+            msprime.PopulationConfiguration(
+                initial_size=sizes[0], metadata=populations[0].asdict()
+            )
+        ],
+    )
+
+
+_species.add_demographic_model(_sma_1pop())
+
+
+def _afr_2epoch():
+    N_A = 746148
+    N_0 = 100218
+    t_1 = 568344
+    populations = [
+        stdpopsim.Population(
+            id="SouthMiddleAtlas",
+            description="Arabidopsis Thaliana South Middle Atlas population",
+        )
+    ]
+    return stdpopsim.DemographicModel(
+        id="African2Epoch_1H18",
+        description="South Middle Atlas African two epoch model",
+        long_description="""
+            Model estimated from site frequency spectrum of synonymous
+            SNPs from African South Middle Atlas samples using
+            Williamson et al. 2005 methodology. Values come from supplementary
+            table 1 of Huber et al 2018. Sizes change from N_A -> N_0 and t_1 is
+            time of the second epoch.
+        """,
+        populations=populations,
+        citations=[
+            stdpopsim.Citation(
+                author="Huber et al.",
+                year=2018,
+                doi="https://doi.org/10.1038/s41467-018-05281-7",
+                reasons={stdpopsim.CiteReason.DEM_MODEL},
+            )
+        ],
+        generation_time=1,
+        population_configurations=[
+            msprime.PopulationConfiguration(
+                initial_size=N_0, metadata=populations[0].asdict()
+            )
+        ],
+        demographic_events=[
+            msprime.PopulationParametersChange(
+                time=t_1, initial_size=N_A, population_id=0
+            )
+        ],
+    )
+
+
+_species.add_demographic_model(_afr_2epoch())
+
+
+def _afr_3epoch():
+    # values from the supplementary Table 1.
+    # the size changed as N_A -> N_2 -> N_3.
+    # t_2 is time of 2nd epoch and t_3 of the third epoch
+    N_A = 161744
+    N_2 = 24076
+    N_3 = 203077
+    t_2 = 7420
+    t_3 = 14534
+    populations = [
+        stdpopsim.Population(
+            id="SouthMiddleAtlas",
+            description="Arabidopsis Thaliana South Middle Atlas population",
+        )
+    ]
+    return stdpopsim.DemographicModel(
+        id="African3Epoch_1H18",
+        description="South Middle Atlas African three epoch model",
+        long_description="""
+            Model estimated from site frequency spectrum of synonymous
+            SNPs from African (South Middle Atlas) samples using Williamson et
+            al. 2005 methodology. Values come from supplementary table 1 of
+            Huber et al 2018. Sizes change from N_A -> N_2 -> N_3 and t_2 is
+            the time of the second epoch and t_3 is the time of the 3rd epoch.
+        """,
+        populations=populations,
+        citations=[
+            stdpopsim.Citation(
+                author="Huber et al.",
+                year=2018,
+                doi="https://doi.org/10.1038/s41467-018-05281-7",
+                reasons={stdpopsim.CiteReason.DEM_MODEL},
+            )
+        ],
+        generation_time=1,
+        population_configurations=[
+            msprime.PopulationConfiguration(
+                initial_size=N_3, metadata=populations[0].asdict()
+            )
+        ],
+        demographic_events=[
+            msprime.PopulationParametersChange(
+                time=t_3, initial_size=N_2, population_id=0
+            ),
+            msprime.PopulationParametersChange(
+                time=t_2 + t_3, initial_size=N_A, population_id=0
+            ),
+        ],
+    )
+
+
+_species.add_demographic_model(_afr_3epoch())
diff --git a/stdpopsim/catalog/AraTha/genetic_maps.py b/stdpopsim/catalog/AraTha/genetic_maps.py
new file mode 100644
index 000000000..0aabed3de
--- /dev/null
+++ b/stdpopsim/catalog/AraTha/genetic_maps.py
@@ -0,0 +1,33 @@
+import stdpopsim
+
+_species = stdpopsim.get_species("AraTha")
+
+
+_gm = stdpopsim.GeneticMap(
+    species=_species,
+    id="SalomeAveraged_TAIR7",
+    description="Crossover frequency map averaged over 17 populations",
+    long_description="""
+        This map is based on the study of crossover frequencies in over 7000
+        plants in 17 F2 populations derived from crosses between 18 A. thaliana
+        accessions. Salomé et al provide genetic maps for each of these
+        populations. To get a single map for each chromosome, the Haldane map
+        function distances were converted to recombination rates (cM/Mb) for
+        each cross and then averaged across the 17 populations using loess.
+        """,
+    url=(
+        "https://stdpopsim.s3-us-west-2.amazonaws.com/genetic_maps/"
+        "AraTha/salome2012_maps.tar.gz"
+    ),
+    sha256="49745e1cab87d59e33eacfdf66303839632d3b07883dd55a99fe1dc27b336ac6",
+    file_pattern="arab_chr{id}_map_loess.txt",
+    citations=[
+        stdpopsim.Citation(
+            doi="https://doi.org/10.1038/hdy.2011.95",
+            author="Salomé et al.",
+            year=2012,
+            reasons={stdpopsim.CiteReason.GEN_MAP},
+        )
+    ],
+)
+_species.add_genetic_map(_gm)
diff --git a/stdpopsim/catalog/AraTha/species.py b/stdpopsim/catalog/AraTha/species.py
new file mode 100644
index 000000000..561ac28dc
--- /dev/null
+++ b/stdpopsim/catalog/AraTha/species.py
@@ -0,0 +1,82 @@
+import stdpopsim
+from . import genome_data
+
+
+_mean_recombination_rate = 200 / 124000 / 2 / 1e6
+
+_recombination_rate_data = {str(j): _mean_recombination_rate for j in range(1, 6)}
+_recombination_rate_data["Mt"] = 0
+_recombination_rate_data["Pt"] = 0  # JK Is this correct??
+
+
+# mutation rate from Ossowski 2010 Science
+# recombination value from Huber et al 2014 MBE
+# rho=200/Mb, assume Ne=124,000, rho=2*Ne*r
+_chromosomes = []
+for name, data in genome_data.data["chromosomes"].items():
+    _chromosomes.append(
+        stdpopsim.Chromosome(
+            id=name,
+            length=data["length"],
+            synonyms=data["synonyms"],
+            mutation_rate=7e-9,
+            recombination_rate=_recombination_rate_data[name],
+        )
+    )
+
+_genome = stdpopsim.Genome(
+    chromosomes=_chromosomes,
+    assembly_name=genome_data.data["assembly_name"],
+    assembly_accession=genome_data.data["assembly_accession"],
+    mutation_rate_citations=[
+        stdpopsim.Citation(
+            author="Ossowski et al.",
+            year="2010",
+            doi="https://doi.org/10.1126/science.1180677",
+            reasons={stdpopsim.CiteReason.MUT_RATE},
+        )
+    ],
+    recombination_rate_citations=[
+        stdpopsim.Citation(
+            author="Huber et al.",
+            year="2014",
+            doi="https://doi.org/10.1093/molbev/msu247",
+            reasons={stdpopsim.CiteReason.REC_RATE},
+        )
+    ],
+    assembly_citations=[
+        stdpopsim.Citation(
+            doi="https://doi.org/10.1093/nar/gkm965",
+            year="2007",
+            author="Swarbreck et al.",
+            reasons={stdpopsim.CiteReason.ASSEMBLY},
+        )
+    ],
+)
+
+_species = stdpopsim.Species(
+    id="AraTha",
+    name="Arabidopsis thaliana",
+    common_name="A. thaliana",
+    genome=_genome,
+    generation_time=1.0,
+    generation_time_citations=[
+        stdpopsim.Citation(
+            doi="https://doi.org/10.1890/0012-9658(2002)083[1006:GTINSO]2.0.CO;2",
+            year="2002",
+            author="Donohue",
+            reasons={stdpopsim.CiteReason.GEN_TIME},
+        )
+    ],
+    population_size=10 ** 4,
+    population_size_citations=[
+        stdpopsim.Citation(
+            doi="https://doi.org/10.1016/j.cell.2016.05.063",
+            year="2016",
+            author="1001GenomesConsortium",
+            reasons={stdpopsim.CiteReason.POP_SIZE},
+        )
+    ],
+)
+
+stdpopsim.register_species(_species)
diff --git a/stdpopsim/catalog/BosTau/__init__.py b/stdpopsim/catalog/BosTau/__init__.py
index de039a4ed..30586d4d5 100644
--- a/stdpopsim/catalog/BosTau/__init__.py
+++ b/stdpopsim/catalog/BosTau/__init__.py
@@ -1,190 +1,5 @@
 """
 Catalog definitions for Bos Taurus
 """
-import collections
-
-import msprime
-
-import stdpopsim
-from . import genome_data
-
-###########################################################
-#
-# Genome definition
-#
-###########################################################
-
-# De novo assembly of the cattle reference genome with single-molecule sequencing.
-_RosenEtAl = stdpopsim.Citation(
-    doi="https://doi.org/10.1093/gigascience/giaa021",
-    year="2020",
-    author="Rosen et al.",
-)
-
-# Frequency of mosaicism points towards mutation-prone early cleavage
-# cell divisions in cattle.
-_HarlandEtAl = stdpopsim.Citation(
-    author="Harland et al.",
-    year="2017",
-    # BioRxiv preprint
-    doi="https://doi.org/10.1101/079863",
-)
-
-# Cattle Sex-Specific Recombination and Genetic Control from a
-# Large Pedigree Analysis.
-_MaEtAl = stdpopsim.Citation(
-    author="Ma et al.",
-    year="2015",
-    doi="https://doi.org/10.1371/journal.pgen.1005387",
-)
-
-# Inferring Demography from Runs of Homozygosity in Whole-Genome Sequence,
-# with Correction for Sequence Errors.
-_MacLeodEtAl = stdpopsim.Citation(
-    doi="https://doi.org/10.1093/molbev/mst125",
-    year="2013",
-    author="MacLeod et al.",
-)
-
-# Recombination rate has been derived from dairy cattle crossovers
-# per meiosis, by taking the average between females and males and then
-# dividing by the whole genome length (equal to the sum of chromosome
-# lengths used above).
-# From Ma et al. (2015), 25.5 crossovers per meiosis in males and
-# 23.2 crossovers per meiosis in females, gives an average of 24.35
-# crossovers per meiosis. The sum of chromosome lengths is 2628394923 bp.
-# 24.35 / 2628394923 = 9.26e-9 per bp per generation.
-_genome_wide_recombination_rate = 9.26e-9
-
-_recombination_rate_data = collections.defaultdict(
-    lambda: _genome_wide_recombination_rate
-)
-# Set some exceptions for non-recombining chrs.
-_recombination_rate_data["MT"] = 0
-
-_chromosomes = []
-for name, data in genome_data.data["chromosomes"].items():
-    _chromosomes.append(
-        stdpopsim.Chromosome(
-            id=name,
-            length=data["length"],
-            synonyms=data["synonyms"],
-            # Harland et al. (2017), sex-averaged estimate per bp per generation.
-            mutation_rate=1.2e-8,
-            recombination_rate=_recombination_rate_data[name],
-        )
-    )
-
-_genome = stdpopsim.Genome(
-    chromosomes=_chromosomes,
-    mutation_rate_citations=[
-        _HarlandEtAl.because(stdpopsim.CiteReason.MUT_RATE),
-    ],
-    recombination_rate_citations=[_MaEtAl.because(stdpopsim.CiteReason.REC_RATE)],
-    assembly_citations=[_RosenEtAl.because(stdpopsim.CiteReason.ASSEMBLY)],
-)
-
-_species = stdpopsim.Species(
-    id="BosTau",
-    name="Bos Taurus",
-    common_name="Cattle",
-    genome=_genome,
-    generation_time=5,
-    generation_time_citations=[_MacLeodEtAl.because(stdpopsim.CiteReason.GEN_TIME)],
-    population_size=62000,
-    population_size_citations=[_MacLeodEtAl.because(stdpopsim.CiteReason.POP_SIZE)],
-)
-
-stdpopsim.register_species(_species)
-
-
-###########################################################
-#
-# Demographic models
-#
-###########################################################
-
-
-def _HolsteinFriesan_1M13():
-    id = "HolsteinFriesian_1M13"
-    description = "Piecewise size changes in Holstein-Friesian cattle."
-    long_description = """
-    The piecewise-constant population size model of Holstein-Friesian cattle
-    from MacLeod et al. 2013. Population sizes were estimated from inferred
-    runs of homozygosity, with parameter values taken from Figure 4A by visual
-    inspection of the plots.
-    """
-    populations = [
-        stdpopsim.Population(id="Holstein_Friesian", description="Holstein Friesian"),
-    ]
-    citations = [_MacLeodEtAl.because(stdpopsim.CiteReason.DEM_MODEL)]
-
-    return stdpopsim.DemographicModel(
-        id=id,
-        description=description,
-        long_description=long_description,
-        populations=populations,
-        citations=citations,
-        generation_time=_species.generation_time,
-        population_configurations=[
-            msprime.PopulationConfiguration(
-                initial_size=90, growth_rate=0.0166, metadata=populations[0].asdict()
-            )
-        ],
-        # Here 'time' should be in generation notation ie. how many
-        # generations ago were that Ne (effective population size)
-        # and growth rate.
-        # Growth rate is "per generation exponential growth rate":
-        # -alpha= [ln(initial_pop_size/next_stage_pop_size)/generation_span_in_years]
-        # For example: ln(90/120)/3= -0.095894024
-        demographic_events=[
-            msprime.PopulationParametersChange(
-                time=1,
-                initial_size=90,
-                growth_rate=-0.095894024,
-                population_id=0,
-            ),  # Ne 90 to 120
-            msprime.PopulationParametersChange(
-                time=4, growth_rate=-0.24465639, population_id=0
-            ),  # Ne 120 to 250
-            msprime.PopulationParametersChange(
-                time=7, growth_rate=-0.0560787, population_id=0
-            ),  # Ne 250 to 350
-            msprime.PopulationParametersChange(
-                time=13, growth_rate=-0.1749704, population_id=0
-            ),  # Ne 350 to 1000
-            msprime.PopulationParametersChange(
-                time=19, growth_rate=-0.0675775, population_id=0
-            ),  # Ne 1000 to 1500
-            msprime.PopulationParametersChange(
-                time=25, growth_rate=-0.0022129, population_id=0
-            ),  # Ne 1500 to 2000
-            msprime.PopulationParametersChange(
-                time=155, growth_rate=-0.0007438, population_id=0
-            ),  # Ne 2000 to 2500
-            msprime.PopulationParametersChange(
-                time=455, growth_rate=-0.0016824, population_id=0
-            ),  # Ne 2500 to 3500
-            msprime.PopulationParametersChange(
-                time=655, growth_rate=-0.0006301, population_id=0
-            ),  # Ne 3500 to 7000
-            msprime.PopulationParametersChange(
-                time=1755, growth_rate=-0.0005945, population_id=0
-            ),  # Ne 7000 to 10000
-            msprime.PopulationParametersChange(
-                time=2355, growth_rate=-0.0005306, population_id=0
-            ),  # Ne 10000 to 17000
-            msprime.PopulationParametersChange(
-                time=3355, growth_rate=-0.0000434, population_id=0
-            ),  # Ne 17000 to 62000
-            msprime.PopulationParametersChange(
-                time=33155, growth_rate=-0.0000, population_id=0
-            ),  # Ne 62000 (model has "coalesced")
-            msprime.PopulationParametersChange(
-                time=933155, growth_rate=-0.0, population_id=0
-            ),
-        ],
-    )
-
-
-_species.add_demographic_model(_HolsteinFriesan_1M13())
+from . import species  # noqa 401
+from . import demographic_models  # noqa 401
diff --git a/stdpopsim/catalog/BosTau/demographic_models.py b/stdpopsim/catalog/BosTau/demographic_models.py
new file mode 100644
index 000000000..e42b761ed
--- /dev/null
+++ b/stdpopsim/catalog/BosTau/demographic_models.py
@@ -0,0 +1,98 @@
+import msprime
+
+import stdpopsim
+
+_species = stdpopsim.get_species("BosTau")
+
+# Inferring Demography from Runs of Homozygosity in Whole-Genome Sequence,
+# with Correction for Sequence Errors.
+_MacLeodEtAl = stdpopsim.Citation(
+    doi="https://doi.org/10.1093/molbev/mst125",
+    year="2013",
+    author="MacLeod et al.",
+)
+
+
+def _HolsteinFriesan_1M13():
+    id = "HolsteinFriesian_1M13"
+    description = "Piecewise size changes in Holstein-Friesian cattle."
+    long_description = """
+    The piecewise-constant population size model of Holstein-Friesian cattle
+    from MacLeod et al. 2013. Population sizes were estimated from inferred
+    runs of homozygosity, with parameter values taken from Figure 4A by visual
+    inspection of the plots.
+    """
+    populations = [
+        stdpopsim.Population(id="Holstein_Friesian", description="Holstein Friesian"),
+    ]
+    citations = [_MacLeodEtAl.because(stdpopsim.CiteReason.DEM_MODEL)]
+
+    return stdpopsim.DemographicModel(
+        id=id,
+        description=description,
+        long_description=long_description,
+        populations=populations,
+        citations=citations,
+        generation_time=_species.generation_time,
+        population_configurations=[
+            msprime.PopulationConfiguration(
+                initial_size=90, growth_rate=0.0166, metadata=populations[0].asdict()
+            )
+        ],
+        # Here 'time' should be in generation notation ie. how many
+        # generations ago were that Ne (effective population size)
+        # and growth rate.
+        # Growth rate is "per generation exponential growth rate":
+        # -alpha= [ln(initial_pop_size/next_stage_pop_size)/generation_span_in_years]
+        # For example: ln(90/120)/3= -0.095894024
+        demographic_events=[
+            msprime.PopulationParametersChange(
+                time=1,
+                initial_size=90,
+                growth_rate=-0.095894024,
+                population_id=0,
+            ),  # Ne 90 to 120
+            msprime.PopulationParametersChange(
+                time=4, growth_rate=-0.24465639, population_id=0
+            ),  # Ne 120 to 250
+            msprime.PopulationParametersChange(
+                time=7, growth_rate=-0.0560787, population_id=0
+            ),  # Ne 250 to 350
+            msprime.PopulationParametersChange(
+                time=13, growth_rate=-0.1749704, population_id=0
+            ),  # Ne 350 to 1000
+            msprime.PopulationParametersChange(
+                time=19, growth_rate=-0.0675775, population_id=0
+            ),  # Ne 1000 to 1500
+            msprime.PopulationParametersChange(
+                time=25, growth_rate=-0.0022129, population_id=0
+            ),  # Ne 1500 to 2000
+            msprime.PopulationParametersChange(
+                time=155, growth_rate=-0.0007438, population_id=0
+            ),  # Ne 2000 to 2500
+            msprime.PopulationParametersChange(
+                time=455, growth_rate=-0.0016824, population_id=0
+            ),  # Ne 2500 to 3500
+            msprime.PopulationParametersChange(
+                time=655, growth_rate=-0.0006301, population_id=0
+            ),  # Ne 3500 to 7000
+            msprime.PopulationParametersChange(
+                time=1755, growth_rate=-0.0005945, population_id=0
+            ),  # Ne 7000 to 10000
+            msprime.PopulationParametersChange(
+                time=2355, growth_rate=-0.0005306, population_id=0
+            ),  # Ne 10000 to 17000
+            msprime.PopulationParametersChange(
+                time=3355, growth_rate=-0.0000434, population_id=0
+            ),  # Ne 17000 to 62000
+            msprime.PopulationParametersChange(
+                time=33155, growth_rate=-0.0000, population_id=0
+            ),  # Ne 62000 (model has "coalesced")
+            msprime.PopulationParametersChange(
+                time=933155, growth_rate=-0.0, population_id=0
+            ),
+        ],
+    )
+
+
+_species.add_demographic_model(_HolsteinFriesan_1M13())
diff --git a/stdpopsim/catalog/BosTau/species.py b/stdpopsim/catalog/BosTau/species.py
new file mode 100644
index 000000000..609de15b5
--- /dev/null
+++ b/stdpopsim/catalog/BosTau/species.py
@@ -0,0 +1,87 @@
+import collections
+
+import stdpopsim
+from . import genome_data
+
+# De novo assembly of the cattle reference genome with single-molecule sequencing.
+_RosenEtAl = stdpopsim.Citation(
+    doi="https://doi.org/10.1093/gigascience/giaa021",
+    year="2020",
+    author="Rosen et al.",
+)
+
+# Frequency of mosaicism points towards mutation-prone early cleavage
+# cell divisions in cattle.
+_HarlandEtAl = stdpopsim.Citation(
+    author="Harland et al.",
+    year="2017",
+    # BioRxiv preprint
+    doi="https://doi.org/10.1101/079863",
+)
+
+# Cattle Sex-Specific Recombination and Genetic Control from a
+# Large Pedigree Analysis.
+_MaEtAl = stdpopsim.Citation(
+    author="Ma et al.",
+    year="2015",
+    doi="https://doi.org/10.1371/journal.pgen.1005387",
+)
+
+# Inferring Demography from Runs of Homozygosity in Whole-Genome Sequence,
+# with Correction for Sequence Errors.
+_MacLeodEtAl = stdpopsim.Citation(
+    doi="https://doi.org/10.1093/molbev/mst125",
+    year="2013",
+    author="MacLeod et al.",
+)
+
+# Recombination rate has been derived from dairy cattle crossovers
+# per meiosis, by taking the average between females and males and then
+# dividing by the whole genome length (equal to the sum of chromosome
+# lengths used above).
+# From Ma et al. (2015), 25.5 crossovers per meiosis in males and
+# 23.2 crossovers per meiosis in females, gives an average of 24.35
+# crossovers per meiosis. The sum of chromosome lengths is 2628394923 bp.
+# 24.35 / 2628394923 = 9.26e-9 per bp per generation.
+_genome_wide_recombination_rate = 9.26e-9
+
+_recombination_rate_data = collections.defaultdict(
+    lambda: _genome_wide_recombination_rate
+)
+# Set some exceptions for non-recombining chrs.
+_recombination_rate_data["MT"] = 0
+
+_chromosomes = []
+for name, data in genome_data.data["chromosomes"].items():
+    _chromosomes.append(
+        stdpopsim.Chromosome(
+            id=name,
+            length=data["length"],
+            synonyms=data["synonyms"],
+            # Harland et al. (2017), sex-averaged estimate per bp per generation.
+            mutation_rate=1.2e-8,
+            recombination_rate=_recombination_rate_data[name],
+        )
+    )
+
+_genome = stdpopsim.Genome(
+    chromosomes=_chromosomes,
+    mutation_rate_citations=[
+        _HarlandEtAl.because(stdpopsim.CiteReason.MUT_RATE),
+    ],
+    recombination_rate_citations=[_MaEtAl.because(stdpopsim.CiteReason.REC_RATE)],
+    assembly_citations=[_RosenEtAl.because(stdpopsim.CiteReason.ASSEMBLY)],
+)
+
+_species = stdpopsim.Species(
+    id="BosTau",
+    name="Bos Taurus",
+    common_name="Cattle",
+    genome=_genome,
+    generation_time=5,
+    generation_time_citations=[_MacLeodEtAl.because(stdpopsim.CiteReason.GEN_TIME)],
+    population_size=62000,
+    population_size_citations=[_MacLeodEtAl.because(stdpopsim.CiteReason.POP_SIZE)],
+)
+
+stdpopsim.register_species(_species)
diff --git a/stdpopsim/catalog/CanFam/__init__.py b/stdpopsim/catalog/CanFam/__init__.py
index 0a3cdcda6..268babc5a 100644
--- a/stdpopsim/catalog/CanFam/__init__.py
+++ b/stdpopsim/catalog/CanFam/__init__.py
@@ -1,148 +1,5 @@
-import stdpopsim
-from . import genome_data
-
-# Recombination rates are the means calculated from the Campbell2016 map.
-# Recombination maps can be found at https://github.com/cflerin/dog_recombination
-# ChrX isn't included in that genetic map, so the weighted genome-wide average is
-# provided here instead.
-_recombination_rate_data = {
-    "1": 7.636001498077e-09,
-    "2": 8.798516284201015e-09,
-    "3": 8.000868549297668e-09,
-    "4": 8.052302321538563e-09,
-    "5": 9.344333839828404e-09,
-    "6": 8.192191165856811e-09,
-    "7": 7.293465322857858e-09,
-    "8": 8.291312822214086e-09,
-    "9": 9.287722798450121e-09,
-    "10": 9.107151930130527e-09,
-    "11": 7.639449804041734e-09,
-    "12": 7.761061128067527e-09,
-    "13": 8.413015225299971e-09,
-    "14": 9.028123091776737e-09,
-    "15": 7.856754982030383e-09,
-    "16": 8.614063697877811e-09,
-    "17": 9.718834385033715e-09,
-    "18": 1.029925446520415e-08,
-    "19": 1.0425051705950442e-08,
-    "20": 9.990971109010329e-09,
-    "21": 1.0339033506095636e-08,
-    "22": 8.61504863805126e-09,
-    "23": 9.1266350068389e-09,
-    "24": 1.1146043069685935e-08,
-    "25": 1.1543746646856006e-08,
-    "26": 1.2084571786110922e-08,
-    "27": 1.1260328390864541e-08,
-    "28": 1.2463587044656355e-08,
-    "29": 1.1013629677730708e-08,
-    "30": 1.1687642441683382e-08,
-    "31": 1.1397713284329192e-08,
-    "32": 1.1555927931648279e-08,
-    "33": 1.3339402745926785e-08,
-    "34": 1.0483812411227089e-08,
-    "35": 1.4299102611645524e-08,
-    "36": 1.187517782077471e-08,
-    "37": 1.3834580623461596e-08,
-    "38": 1.4363726512881696e-08,
-    "X": 9.506483722244087e-09,
-    "MT": 0,
-}
-
-_LindbladTohEtAl = stdpopsim.Citation(
-    # Genome sequence, comparative analysis and haplotype structure of the
-    # domestic dog.
-    author="Lindblad-Toh et al.",
-    year=2005,
-    doi="https://doi.org/10.1038/nature04338",
-)
-
-_SkoglundEtAl = stdpopsim.Citation(
-    # Ancient wolf genome reveals an early divergence of domestic dog
-    # ancestors and admixture into high-latitude breeds.
-    author="Skoglund et al.",
-    year=2015,
-    doi="https://doi.org/10.1016/j.cub.2015.04.019",
-)
-
-_FranzEtAl = stdpopsim.Citation(
-    # Genomic and archaeological evidence suggest a dual origin of
-    # domestic dogs.
-    author="Franz et al.",
-    year=2016,
-    doi="https://doi.org/10.1126/science.aaf3161",
-)
-
-_CampbellEtAl = stdpopsim.Citation(
-    # A Pedigree-Based Map of Recombination in the Domestic Dog Genome.
-    author="Campbell et al.",
-    year=2016,
-    doi="https://doi.org/10.1534/g3.116.034678",
-)
-
-_chromosomes = []
-for name, data in genome_data.data["chromosomes"].items():
-    _chromosomes.append(
-        stdpopsim.Chromosome(
-            id=name,
-            length=data["length"],
-            synonyms=data["synonyms"],
-            mutation_rate=4e-9,  # based on non-CpG sites only
-            recombination_rate=_recombination_rate_data[name],
-        )
-    )
-
-_genome = stdpopsim.Genome(
-    chromosomes=_chromosomes,
-    assembly_name=genome_data.data["assembly_name"],
-    assembly_accession=genome_data.data["assembly_accession"],
-    mutation_rate_citations=[
-        _SkoglundEtAl.because(stdpopsim.CiteReason.MUT_RATE),
-        _FranzEtAl.because(stdpopsim.CiteReason.MUT_RATE),
-    ],
-    recombination_rate_citations=[_CampbellEtAl.because(stdpopsim.CiteReason.REC_RATE)],
-    assembly_citations=[_LindbladTohEtAl.because(stdpopsim.CiteReason.ASSEMBLY)],
-)
-
-_species = stdpopsim.Species(
-    id="CanFam",
-    name="Canis familiaris",
-    common_name="Dog",
-    genome=_genome,
-    generation_time=3,
-    generation_time_citations=[
-        # Everyone uses 3 years because everyone else uses it.
-        # It's likely higher, at least in wolves:
-        # https://academic.oup.com/mbe/article/35/6/1366/4990884
-        # Reasoning behind a generation time of 3 years:
-        # Consider two use cases for CanFam simulations:
-        # (1) for domestic dog simulations, and (2) for wolf+dog simulations
-        # (or ancestral dogs).
-        # In case (1), maybe 3 year generations are more appropriate because of human
-        # intervention in breeding. In case (2), you might want to match what other
-        # studies have done (thus using 3 year generations), or you might want to
-        # consider what is known about modern wolves.
-    ],
-    population_size=13000,  # ancestral dog size
-    population_size_citations=[_LindbladTohEtAl.because(stdpopsim.CiteReason.POP_SIZE)],
-)
-
-stdpopsim.register_species(_species)
-
-_gm = stdpopsim.GeneticMap(
-    species=_species,
-    id="Campbell2016_CanFam3_1",
-    description="Pedigree-based crossover map from 237 individuals",
-    long_description="""
-        Sex-averaged crossover frequency map based on 163,400 autosomal SNPs
-        genotyped in a pedigree of 237 Labrador Retriever x Greyhound crosses.
-        Genotypes were phased without respect to the pedigree, using SHAPEIT2,
-        recombinations were called using duoHMM, and genetic distances were
-        obtained using Haldane's map function.
-        """,
-    url="https://stdpopsim.s3-us-west-2.amazonaws.com/genetic_maps/"
-    "CanFam/dog_genetic_maps.tar.gz",
-    sha256="585afb424615e2fb0825d807db0b10fe1c797a6dbb804ecbb3fef5e8387d194f",
-    file_pattern="chr{id}_average_canFam3.1.txt",
-    citations=[_CampbellEtAl.because(stdpopsim.CiteReason.GEN_MAP)],
-)
-_species.add_genetic_map(_gm)
+"""
+Catalog data for Canis familiaris.
+"""
+from . import species  # noqa 401
+from . import genetic_maps  # noqa 401
diff --git a/stdpopsim/catalog/CanFam/genetic_maps.py b/stdpopsim/catalog/CanFam/genetic_maps.py
new file mode 100644
index 000000000..393ae3835
--- /dev/null
+++ b/stdpopsim/catalog/CanFam/genetic_maps.py
@@ -0,0 +1,29 @@
+import stdpopsim
+
+_species = stdpopsim.get_species("CanFam")
+
+_CampbellEtAl = stdpopsim.Citation(
+    # A Pedigree-Based Map of Recombination in the Domestic Dog Genome.
+    author="Campbell et al.",
+    year=2016,
+    doi="https://doi.org/10.1534/g3.116.034678",
+)
+
+_gm = stdpopsim.GeneticMap(
+    species=_species,
+    id="Campbell2016_CanFam3_1",
+    description="Pedigree-based crossover map from 237 individuals",
+    long_description="""
+        Sex-averaged crossover frequency map based on 163,400 autosomal SNPs
+        genotyped in a pedigree of 237 Labrador Retriever x Greyhound crosses.
+        Genotypes were phased without respect to the pedigree, using SHAPEIT2,
+        recombinations were called using duoHMM, and genetic distances were
+        obtained using Haldane's map function.
+        """,
+    url="https://stdpopsim.s3-us-west-2.amazonaws.com/genetic_maps/"
+    "CanFam/dog_genetic_maps.tar.gz",
+    sha256="585afb424615e2fb0825d807db0b10fe1c797a6dbb804ecbb3fef5e8387d194f",
+    file_pattern="chr{id}_average_canFam3.1.txt",
+    citations=[_CampbellEtAl.because(stdpopsim.CiteReason.GEN_MAP)],
+)
+_species.add_genetic_map(_gm)
diff --git a/stdpopsim/catalog/CanFam/species.py b/stdpopsim/catalog/CanFam/species.py
new file mode 100644
index 000000000..ddea23514
--- /dev/null
+++ b/stdpopsim/catalog/CanFam/species.py
@@ -0,0 +1,129 @@
+import stdpopsim
+from . import genome_data
+
+# Recombination rates are the means calculated from the Campbell2016 map.
+# Recombination maps can be found at https://github.com/cflerin/dog_recombination
+# ChrX isn't included in that genetic map, so the weighted genome-wide average is
+# provided here instead.
+_recombination_rate_data = {
+    "1": 7.636001498077e-09,
+    "2": 8.798516284201015e-09,
+    "3": 8.000868549297668e-09,
+    "4": 8.052302321538563e-09,
+    "5": 9.344333839828404e-09,
+    "6": 8.192191165856811e-09,
+    "7": 7.293465322857858e-09,
+    "8": 8.291312822214086e-09,
+    "9": 9.287722798450121e-09,
+    "10": 9.107151930130527e-09,
+    "11": 7.639449804041734e-09,
+    "12": 7.761061128067527e-09,
+    "13": 8.413015225299971e-09,
+    "14": 9.028123091776737e-09,
+    "15": 7.856754982030383e-09,
+    "16": 8.614063697877811e-09,
+    "17": 9.718834385033715e-09,
+    "18": 1.029925446520415e-08,
+    "19": 1.0425051705950442e-08,
+    "20": 9.990971109010329e-09,
+    "21": 1.0339033506095636e-08,
+    "22": 8.61504863805126e-09,
+    "23": 9.1266350068389e-09,
+    "24": 1.1146043069685935e-08,
+    "25": 1.1543746646856006e-08,
+    "26": 1.2084571786110922e-08,
+    "27": 1.1260328390864541e-08,
+    "28": 1.2463587044656355e-08,
+    "29": 1.1013629677730708e-08,
+    "30": 1.1687642441683382e-08,
+    "31": 1.1397713284329192e-08,
+    "32": 1.1555927931648279e-08,
+    "33": 1.3339402745926785e-08,
+    "34": 1.0483812411227089e-08,
+    "35": 1.4299102611645524e-08,
+    "36": 1.187517782077471e-08,
+    "37": 1.3834580623461596e-08,
+    "38": 1.4363726512881696e-08,
+    "X": 9.506483722244087e-09,
+    "MT": 0,
+}
+
+_LindbladTohEtAl = stdpopsim.Citation(
+    # Genome sequence, comparative analysis and haplotype structure of the
+    # domestic dog.
+    author="Lindblad-Toh et al.",
+    year=2005,
+    doi="https://doi.org/10.1038/nature04338",
+)
+
+_SkoglundEtAl = stdpopsim.Citation(
+    # Ancient wolf genome reveals an early divergence of domestic dog
+    # ancestors and admixture into high-latitude breeds.
+    author="Skoglund et al.",
+    year=2015,
+    doi="https://doi.org/10.1016/j.cub.2015.04.019",
+)
+
+_FranzEtAl = stdpopsim.Citation(
+    # Genomic and archaeological evidence suggest a dual origin of
+    # domestic dogs.
+    author="Franz et al.",
+    year=2016,
+    doi="https://doi.org/10.1126/science.aaf3161",
+)
+
+_CampbellEtAl = stdpopsim.Citation(
+    # A Pedigree-Based Map of Recombination in the Domestic Dog Genome.
+    author="Campbell et al.",
+    year=2016,
+    doi="https://doi.org/10.1534/g3.116.034678",
+)
+
+_chromosomes = []
+for name, data in genome_data.data["chromosomes"].items():
+    _chromosomes.append(
+        stdpopsim.Chromosome(
+            id=name,
+            length=data["length"],
+            synonyms=data["synonyms"],
+            mutation_rate=4e-9,  # based on non-CpG sites only
+            recombination_rate=_recombination_rate_data[name],
+        )
+    )
+
+_genome = stdpopsim.Genome(
+    chromosomes=_chromosomes,
+    assembly_name=genome_data.data["assembly_name"],
+    assembly_accession=genome_data.data["assembly_accession"],
+    mutation_rate_citations=[
+        _SkoglundEtAl.because(stdpopsim.CiteReason.MUT_RATE),
+        _FranzEtAl.because(stdpopsim.CiteReason.MUT_RATE),
+    ],
+    recombination_rate_citations=[_CampbellEtAl.because(stdpopsim.CiteReason.REC_RATE)],
+    assembly_citations=[_LindbladTohEtAl.because(stdpopsim.CiteReason.ASSEMBLY)],
+)
+
+_species = stdpopsim.Species(
+    id="CanFam",
+    name="Canis familiaris",
+    common_name="Dog",
+    genome=_genome,
+    generation_time=3,
+    generation_time_citations=[
+        # Everyone uses 3 years because everyone else uses it.
+        # It's likely higher, at least in wolves:
+        # https://academic.oup.com/mbe/article/35/6/1366/4990884
+        # Reasoning behind a generation time of 3 years:
+        # Consider two use cases for CanFam simulations:
+        # (1) for domestic dog simulations, and (2) for wolf+dog simulations
+        # (or ancestral dogs).
+        # In case (1), maybe 3 year generations are more appropriate because of human
+        # intervention in breeding. In case (2), you might want to match what other
+        # studies have done (thus using 3 year generations), or you might want to
+        # consider what is known about modern wolves.
+    ],
+    population_size=13000,  # ancestral dog size
+    population_size_citations=[_LindbladTohEtAl.because(stdpopsim.CiteReason.POP_SIZE)],
+)
+
+stdpopsim.register_species(_species)
diff --git a/stdpopsim/catalog/DroMel/__init__.py b/stdpopsim/catalog/DroMel/__init__.py
index 6cfd10782..8ed4e71c3 100644
--- a/stdpopsim/catalog/DroMel/__init__.py
+++ b/stdpopsim/catalog/DroMel/__init__.py
@@ -1,287 +1,6 @@
 """
-Genome, genetic map and demographic model definitions for humans.
+Catalog data for Drosophila Melanogaster
 """
-import collections
-
-import msprime
-
-import stdpopsim
-from . import genome_data
-
-###########################################################
-#
-# Genome definition
-#
-###########################################################
-
-# citations
-_LiAndStephan = stdpopsim.Citation(
-    author="Li et al.", year=2006, doi="https://doi.org/10.1371/journal.pgen.0020166"
-)
-
-_SchriderEtAl = stdpopsim.Citation(
-    author="Schrider et al.",
-    year=2013,
-    doi="https://doi.org/10.1534/genetics.113.151670",
-)
-
-_DosSantosEtAl = stdpopsim.Citation(
-    doi="https://doi.org/10.1093/nar/gku1099",
-    year="2015",
-    author="dos Santos et al.",
-    reasons={stdpopsim.CiteReason.ASSEMBLY},
-)
-
-_genome_wide_estimate = 8.4e-9  # WRONG, underestimate used in S&S!
-
-_recombination_rate_data = collections.defaultdict(lambda: _genome_wide_estimate)
-# Set some exceptions for non-recombining chrs.
-_recombination_rate_data["Y"] = 0
-_recombination_rate_data["mitochondrion_genome"] = 0
-
-_chromosomes = []
-for name, data in genome_data.data["chromosomes"].items():
-    _chromosomes.append(
-        stdpopsim.Chromosome(
-            id=name,
-            length=data["length"],
-            synonyms=data["synonyms"],
-            mutation_rate=5.49e-9,  # citation: _SchriderEtAl
-            recombination_rate=_recombination_rate_data[name],
-        )
-    )
-
-
-_genome = stdpopsim.Genome(
-    chromosomes=_chromosomes,
-    assembly_name=genome_data.data["assembly_name"],
-    assembly_accession=genome_data.data["assembly_accession"],
-    mutation_rate_citations=[_SchriderEtAl.because(stdpopsim.CiteReason.MUT_RATE)],
-    assembly_citations=[_DosSantosEtAl],
-)
-
-_species = stdpopsim.Species(
-    id="DroMel",
-    name="Drosophila melanogaster",
-    common_name="D. melanogaster",
-    genome=_genome,
-    generation_time=0.1,
-    generation_time_citations=[_LiAndStephan.because(stdpopsim.CiteReason.GEN_TIME)],
-    population_size=1720600,
-    population_size_citations=[_LiAndStephan.because(stdpopsim.CiteReason.POP_SIZE)],
-)
-
-stdpopsim.register_species(_species)
-
-
-###########################################################
-#
-# Genetic maps
-#
-###########################################################
-
-
-_gm = stdpopsim.GeneticMap(
-    species=_species,
-    id="ComeronCrossover_dm6",
-    description="Crossover map from meioses products of 8 lab crosses",
-    long_description="""
-        The crossover map from a study of 8 crosses of 12 highly
-        inbred lines of D. melanogaster. This is based on the
-        products of 5,860 female meioses from whole genome sequencing data.
-        Recombination rates were calculated from the density of individual
-        recombination events that were detected in crosses. This map was
-        subsequently lifted over to the dm6 assembly.
-        """,
-    url=(
-        "https://stdpopsim.s3-us-west-2.amazonaws.com/genetic_maps/"
-        "DroMel/comeron2012_maps.tar.gz"
-    ),
-    sha256="08185a0e3b0ad26eefe69fc6bdb8f3f599a760e11e87dd343335b33d1563f62a",
-    file_pattern="genetic_map_comeron2012_dm6_chr{id}.txt",
-    citations=[
-        stdpopsim.Citation(
-            author="Comeron et al",
-            doi="https://doi.org/10.1371/journal.pgen.1002905",
-            year=2012,
-            reasons={stdpopsim.CiteReason.GEN_MAP},
-        )
-    ],
-)
-_species.add_genetic_map(_gm)
-
-_gm = stdpopsim.GeneticMap(
-    species=_species,
-    id="ComeronCrossoverV2_dm6",
-    description="Crossover map from meioses products of 8 lab crosses",
-    long_description="""
-        The crossover map from a study of 8 crosses of 12 highly
-        inbred lines of D. melanogaster. This is based on the
-        products of 5,860 female meioses from whole genome sequencing data.
-        Recombination rates were calculated from the density of individual
-        recombination events that were detected in crosses. This map was
-        subsequently lifted over to the dm6 assembly using the available
-        maintenance code command:
-
-        python liftOver_comeron2012.py \
-            --winLen 1000 \
-            --gapThresh 1000000 \
-            --useAdjacentAvg \
-            --retainIntermediates
-        """,
-    url=(
-        "https://stdpopsim.s3-us-west-2.amazonaws.com/genetic_maps/"
-        "DroMel/comeron2012v2_maps.tar.gz"
-    ),
-    sha256="79484e6042c36b06946af672cabf2ec3012b6bfc6e0018517c8e7be4769cd334",
-    file_pattern="genetic_map_comeron2012v2_dm6_chr{id}.txt",
-    citations=[
-        stdpopsim.Citation(
-            author="Comeron et al",
-            doi="https://doi.org/10.1371/journal.pgen.1002905",
-            year=2012,
-            reasons={stdpopsim.CiteReason.GEN_MAP},
-        )
-    ],
-)
-_species.add_genetic_map(_gm)
-
-
-###########################################################
-#
-# Demographic models
-#
-###########################################################
-
-
-# population definitions that are reused.
-_afr_population = stdpopsim.Population(
-    id="AFR", description="African D. melanogaster population"
-)
-_eur_population = stdpopsim.Population(
-    id="EUR", description="European D. melanogaster population"
-)
-
-
-def _afr_3epoch():
-    id = "African3Epoch_1S16"
-    description = "Three epoch African population"
-    long_description = """
-        The three epoch (modern, bottleneck, ancestral) model estimated for a
-        single African Drosophila Melanogaster population from Sheehan and Song (2016).
-        Population sizes are estimated by a
-        deep learning model trained on simulation data. NOTE: Due to differences in
-        coalescence units between PSMC (2N) and msms (4N) the number of generations were
-        doubled from PSMC estimates when simulating data from msms in the original
-        publication. We have faithfully represented the published model here.
-    """
-    populations = [_afr_population]
-    citations = [
-        stdpopsim.Citation(
-            author="Sheehan and Song",
-            year=2016,
-            doi="https://doi.org/10.1371/journal.pcbi.1004845",
-            reasons={stdpopsim.CiteReason.DEM_MODEL},
-        )
-    ]
-    generation_time = _species.generation_time
-    citations.extend(_species.generation_time_citations)
-
-    # Parameter values from "Simulating Data" section
-    # these are assumptions, not estimates
-    N_ref = 100000
-    t_1_coal = 0.5
-    t_2_coal = 5.0
-    # estimates from the ANN
-    N_R = 544200
-    N_B = 145300
-    N_A = 652700
-    # Times are provided in 4N_ref generations, so we convert into generations.
-    # generation_time = 10 / year
-    t_1 = t_1_coal * 4 * N_ref
-    t_2 = (t_1_coal + t_2_coal) * 4 * N_ref
-
-    return stdpopsim.DemographicModel(
-        id=id,
-        description=description,
-        long_description=long_description,
-        populations=populations,
-        citations=citations,
-        generation_time=generation_time,
-        population_configurations=[
-            msprime.PopulationConfiguration(
-                initial_size=N_R, metadata=populations[0].asdict()
-            )
-        ],
-        demographic_events=[
-            # Size change at bottleneck (back in time; BIT)
-            msprime.PopulationParametersChange(
-                time=t_1, initial_size=N_B, population_id=0
-            ),
-            # Size change at recovery (BIT)
-            msprime.PopulationParametersChange(
-                time=t_2, initial_size=N_A, population_id=0
-            ),
-        ],
-    )
-
-
-_species.add_demographic_model(_afr_3epoch())
-
-
-def _ooa_2():
-    id = "OutOfAfrica_2L06"
-    description = "Three epoch model for African and European populations"
-    long_description = """
-        The three epoch (modern, bottleneck, ancestral) model estimated for two
-        Drosophila Melanogaster populations: African (ancestral) and European (derived)
-        from Li and Stephan (2006).
-    """
-    populations = [_afr_population, _eur_population]
-    citations = [_LiAndStephan.because(stdpopsim.CiteReason.DEM_MODEL)]
-    generation_time = _species.generation_time
-    citations.extend(_species.generation_time_citations)
-
-    # African Parameter values from "Demographic History of the African
-    # Population" section
-    N_A0 = 8.603e06
-    t_A0 = 600000  # assuming 10 generations / year
-    N_A1 = N_A0 / 5.0
-
-    # European Parameter values from "Demo History of Euro Population"
-    N_E0 = 1.075e06
-    N_E1 = 2200
-    t_AE = 158000  # generations
-    t_E1 = t_AE - 3400
-
-    return stdpopsim.DemographicModel(
-        id=id,
-        description=description,
-        long_description=long_description,
-        populations=populations,
-        citations=citations,
-        generation_time=generation_time,
-        population_configurations=[
-            msprime.PopulationConfiguration(
-                initial_size=N_A0, metadata=populations[0].asdict()
-            ),
-            msprime.PopulationConfiguration(
-                initial_size=N_E0, metadata=populations[1].asdict()
-            ),
-        ],
-        demographic_events=[
-            # Size change at Euro bottleneck
-            msprime.PopulationParametersChange(
-                time=t_E1, initial_size=N_E1, population_id=1
-            ),
-            # Split
-            msprime.MassMigration(time=t_AE, source=1, destination=0, proportion=1.0),
-            # African bottleneck
-            msprime.PopulationParametersChange(
-                time=t_A0, initial_size=N_A1, population_id=0
-            ),
-        ],
-    )
-
-
-_species.add_demographic_model(_ooa_2())
+from . import species  # noqa 401
+from . import genetic_maps  # noqa 401
+from . import demographic_models  # noqa 401
diff --git a/stdpopsim/catalog/DroMel/demographic_models.py b/stdpopsim/catalog/DroMel/demographic_models.py
new file mode 100644
index 000000000..27c2aa9f2
--- /dev/null
+++ b/stdpopsim/catalog/DroMel/demographic_models.py
@@ -0,0 +1,142 @@
+import msprime
+
+import stdpopsim
+
+_species = stdpopsim.get_species("DroMel")
+
+
+_LiAndStephan = stdpopsim.Citation(
+    author="Li et al.", year=2006, doi="https://doi.org/10.1371/journal.pgen.0020166"
+)
+
+# population definitions that are reused.
+_afr_population = stdpopsim.Population(
+    id="AFR", description="African D. melanogaster population"
+)
+_eur_population = stdpopsim.Population(
+    id="EUR", description="European D. melanogaster population"
+)
+
+
+def _afr_3epoch():
+    id = "African3Epoch_1S16"
+    description = "Three epoch African population"
+    long_description = """
+        The three epoch (modern, bottleneck, ancestral) model estimated for a
+        single African Drosophila Melanogaster population from Sheehan and Song (2016).
+        Population sizes are estimated by a
+        deep learning model trained on simulation data. NOTE: Due to differences in
+        coalescence units between PSMC (2N) and msms (4N) the number of generations were
+        doubled from PSMC estimates when simulating data from msms in the original
+        publication. We have faithfully represented the published model here.
+    """
+    populations = [_afr_population]
+    citations = [
+        stdpopsim.Citation(
+            author="Sheehan and Song",
+            year=2016,
+            doi="https://doi.org/10.1371/journal.pcbi.1004845",
+            reasons={stdpopsim.CiteReason.DEM_MODEL},
+        )
+    ]
+    generation_time = _species.generation_time
+    citations.extend(_species.generation_time_citations)
+
+    # Parameter values from "Simulating Data" section
+    # these are assumptions, not estimates
+    N_ref = 100000
+    t_1_coal = 0.5
+    t_2_coal = 5.0
+    # estimates from the ANN
+    N_R = 544200
+    N_B = 145300
+    N_A = 652700
+    # Times are provided in 4N_ref generations, so we convert into generations.
+    # generation_time = 10 / year
+    t_1 = t_1_coal * 4 * N_ref
+    t_2 = (t_1_coal + t_2_coal) * 4 * N_ref
+
+    return stdpopsim.DemographicModel(
+        id=id,
+        description=description,
+        long_description=long_description,
+        populations=populations,
+        citations=citations,
+        generation_time=generation_time,
+        population_configurations=[
+            msprime.PopulationConfiguration(
+                initial_size=N_R, metadata=populations[0].asdict()
+            )
+        ],
+        demographic_events=[
+            # Size change at bottleneck (back in time; BIT)
+            msprime.PopulationParametersChange(
+                time=t_1, initial_size=N_B, population_id=0
+            ),
+            # Size change at recovery (BIT)
+            msprime.PopulationParametersChange(
+                time=t_2, initial_size=N_A, population_id=0
+            ),
+        ],
+    )
+
+
+_species.add_demographic_model(_afr_3epoch())
+
+
+def _ooa_2():
+    id = "OutOfAfrica_2L06"
+    description = "Three epoch model for African and European populations"
+    long_description = """
+        The three epoch (modern, bottleneck, ancestral) model estimated for two
+        Drosophila Melanogaster populations: African (ancestral) and European (derived)
+        from Li and Stephan (2006).
+    """
+    populations = [_afr_population, _eur_population]
+    citations = [_LiAndStephan.because(stdpopsim.CiteReason.DEM_MODEL)]
+    generation_time = _species.generation_time
+    citations.extend(_species.generation_time_citations)
+
+    # African Parameter values from "Demographic History of the African
+    # Population" section
+    N_A0 = 8.603e06
+    t_A0 = 600000  # assuming 10 generations / year
+    N_A1 = N_A0 / 5.0
+
+    # European Parameter values from "Demo History of Euro Population"
+    N_E0 = 1.075e06
+    N_E1 = 2200
+    t_AE = 158000  # generations
+    t_E1 = t_AE - 3400
+
+    return stdpopsim.DemographicModel(
+        id=id,
+        description=description,
+        long_description=long_description,
+        populations=populations,
+        citations=citations,
+        generation_time=generation_time,
+        population_configurations=[
+            msprime.PopulationConfiguration(
+                initial_size=N_A0, metadata=populations[0].asdict()
+            ),
+            msprime.PopulationConfiguration(
+                initial_size=N_E0, metadata=populations[1].asdict()
+            ),
+        ],
+        demographic_events=[
+            # Size change at Euro bottleneck
+            msprime.PopulationParametersChange(
+                time=t_E1, initial_size=N_E1, population_id=1
+            ),
+            # Split
+            msprime.MassMigration(time=t_AE, source=1, destination=0, proportion=1.0),
+            # African bottleneck
+            msprime.PopulationParametersChange(
+                time=t_A0, initial_size=N_A1, population_id=0
+            ),
+        ],
+    )
+
+
+_species.add_demographic_model(_ooa_2())
diff --git a/stdpopsim/catalog/DroMel/genetic_maps.py b/stdpopsim/catalog/DroMel/genetic_maps.py
new file mode 100644
index 000000000..de515d9dc
--- /dev/null
+++ b/stdpopsim/catalog/DroMel/genetic_maps.py
@@ -0,0 +1,68 @@
+import stdpopsim
+
+_species = stdpopsim.get_species("DroMel")
+
+_gm = stdpopsim.GeneticMap(
+    species=_species,
+    id="ComeronCrossover_dm6",
+    description="Crossover map from meioses products of 8 lab crosses",
+    long_description="""
+        The crossover map from a study of 8 crosses of 12 highly
+        inbred lines of D. melanogaster. This is based on the
+        products of 5,860 female meioses from whole genome sequencing data.
+        Recombination rates were calculated from the density of individual
+        recombination events that were detected in crosses. This map was
+        subsequently lifted over to the dm6 assembly.
+        """,
+    url=(
+        "https://stdpopsim.s3-us-west-2.amazonaws.com/genetic_maps/"
+        "DroMel/comeron2012_maps.tar.gz"
+    ),
+    sha256="08185a0e3b0ad26eefe69fc6bdb8f3f599a760e11e87dd343335b33d1563f62a",
+    file_pattern="genetic_map_comeron2012_dm6_chr{id}.txt",
+    citations=[
+        stdpopsim.Citation(
+            author="Comeron et al",
+            doi="https://doi.org/10.1371/journal.pgen.1002905",
+            year=2012,
+            reasons={stdpopsim.CiteReason.GEN_MAP},
+        )
+    ],
+)
+_species.add_genetic_map(_gm)
+
+_gm = stdpopsim.GeneticMap(
+    species=_species,
+    id="ComeronCrossoverV2_dm6",
+    description="Crossover map from meioses products of 8 lab crosses",
+    long_description="""
+        The crossover map from a study of 8 crosses of 12 highly
+        inbred lines of D. melanogaster. This is based on the
+        products of 5,860 female meioses from whole genome sequencing data.
+        Recombination rates were calculated from the density of individual
+        recombination events that were detected in crosses. This map was
+        subsequently lifted over to the dm6 assembly using the available
+        maintenance code command:
+
+        python liftOver_comeron2012.py \
+            --winLen 1000 \
+            --gapThresh 1000000 \
+            --useAdjacentAvg \
+            --retainIntermediates
+        """,
+    url=(
+        "https://stdpopsim.s3-us-west-2.amazonaws.com/genetic_maps/"
+        "DroMel/comeron2012v2_maps.tar.gz"
+    ),
+    sha256="79484e6042c36b06946af672cabf2ec3012b6bfc6e0018517c8e7be4769cd334",
+    file_pattern="genetic_map_comeron2012v2_dm6_chr{id}.txt",
+    citations=[
+        stdpopsim.Citation(
+            author="Comeron et al",
+            doi="https://doi.org/10.1371/journal.pgen.1002905",
+            year=2012,
+            reasons={stdpopsim.CiteReason.GEN_MAP},
+        )
+    ],
+)
+_species.add_genetic_map(_gm)
diff --git a/stdpopsim/catalog/DroMel/species.py b/stdpopsim/catalog/DroMel/species.py
new file mode 100644
index 000000000..ecc501397
--- /dev/null
+++ b/stdpopsim/catalog/DroMel/species.py
@@ -0,0 +1,68 @@
+import collections
+
+import stdpopsim
+from . import genome_data
+
+###########################################################
+#
+# Genome definition
+#
+###########################################################
+
+_LiAndStephan = stdpopsim.Citation(
+    author="Li et al.", year=2006, doi="https://doi.org/10.1371/journal.pgen.0020166"
+)
+
+_SchriderEtAl = stdpopsim.Citation(
+    author="Schrider et al.",
+    year=2013,
+    doi="https://doi.org/10.1534/genetics.113.151670",
+)
+
+_DosSantosEtAl = stdpopsim.Citation(
+    doi="https://doi.org/10.1093/nar/gku1099",
+    year="2015",
+    author="dos Santos et al.",
+    reasons={stdpopsim.CiteReason.ASSEMBLY},
+)
+
+_genome_wide_estimate = 8.4e-9  # WRONG, underestimate used in S&S!
+
+_recombination_rate_data = collections.defaultdict(lambda: _genome_wide_estimate)
+# Set some exceptions for non-recombining chrs.
+_recombination_rate_data["Y"] = 0
+_recombination_rate_data["mitochondrion_genome"] = 0
+
+_chromosomes = []
+for name, data in genome_data.data["chromosomes"].items():
+    _chromosomes.append(
+        stdpopsim.Chromosome(
+            id=name,
+            length=data["length"],
+            synonyms=data["synonyms"],
+            mutation_rate=5.49e-9,  # citation: _SchriderEtAl
+            recombination_rate=_recombination_rate_data[name],
+        )
+    )
+
+
+_genome = stdpopsim.Genome(
+    chromosomes=_chromosomes,
+    assembly_name=genome_data.data["assembly_name"],
+    assembly_accession=genome_data.data["assembly_accession"],
+    mutation_rate_citations=[_SchriderEtAl.because(stdpopsim.CiteReason.MUT_RATE)],
+    assembly_citations=[_DosSantosEtAl],
+)
+
+_species = stdpopsim.Species(
+    id="DroMel",
+    name="Drosophila melanogaster",
+    common_name="D. melanogaster",
+    genome=_genome,
+    generation_time=0.1,
+    generation_time_citations=[_LiAndStephan.because(stdpopsim.CiteReason.GEN_TIME)],
+    population_size=1720600,
+    population_size_citations=[_LiAndStephan.because(stdpopsim.CiteReason.POP_SIZE)],
+)
+
+stdpopsim.register_species(_species)
diff --git a/stdpopsim/catalog/EscCol/__init__.py b/stdpopsim/catalog/EscCol/__init__.py
index 2931bccd3..8266c1ffa 100644
--- a/stdpopsim/catalog/EscCol/__init__.py
+++ b/stdpopsim/catalog/EscCol/__init__.py
@@ -1,79 +1,4 @@
 """
 Genome and demographic model definitions for Escherichia coli.
 """
-import stdpopsim
-from . import genome_data
-
-###########################################################
-#
-# Genome definition
-#
-###########################################################
-
-_hartl_et_al = stdpopsim.Citation(
-    author="Hartl, Moriyama, and Sawyer",
-    year="1994",
-    # doesn't have a doi
-    doi="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1206133/",
-)
-
-_sezonov_et_al = stdpopsim.Citation(
-    author="Sezonov et al.", year="2007", doi="https://doi.org/10.1128/JB.01368-07"
-)
-
-_wielgoss_et_al = stdpopsim.Citation(
-    author="Wielgoss et al.", year="2011", doi="https://doi.org/10.1534/g3.111.000406"
-)
-
-_blattner_et_al = stdpopsim.Citation(
-    author="Blattner et al.", year="1997", doi="10.1126/science.277.5331.1453"
-)
-
-_chromosomes = []
-for name, data in genome_data.data["chromosomes"].items():
-    _chromosomes.append(
-        stdpopsim.Chromosome(
-            id=name,
-            length=data["length"],
-            synonyms=data["synonyms"],
-            # Wielgoss et al. (2011) calculated for strain REL606,
-            # from synonymous substitutions over 40,000 generations.
-            mutation_rate=8.9e-11,
-            recombination_rate=0.0,
-        )
-    )
-
-# mean_conversion_rate=8.9e-11 # not implemented yet!
-# mean_conversion_length=542 # not implemented yet!
-
-_genome = stdpopsim.Genome(
-    chromosomes=_chromosomes,
-    assembly_name=genome_data.data["assembly_name"],
-    assembly_accession=genome_data.data["assembly_accession"],
-    mutation_rate_citations=[
-        _wielgoss_et_al.because(stdpopsim.CiteReason.MUT_RATE),
-    ],
-    assembly_citations=[_blattner_et_al.because(stdpopsim.CiteReason.ASSEMBLY)],
-)
-
-
-_species = stdpopsim.Species(
-    id="EscCol",
-    name="Escherichia coli",
-    common_name="E. coli",
-    # We use the K-12 strain, because the parameters we're using more
-    # closely match this strain than the ensembl default (HUSEC2011).
-    ensembl_id="escherichia_coli_str_k_12_substr_mg1655",
-    genome=_genome,
-    # E. coli K-12 strain MG1655 "doubling time during steady-state growth in
-    # Luria-Bertani broth was 20 min".
-    generation_time=0.00003805175,  # 1.0 / (525600 min/year / 20 min/gen)
-    generation_time_citations=[_sezonov_et_al.because(stdpopsim.CiteReason.GEN_TIME)],
-    # Hartl et al. calculated Ne for "natural isolates of E. coli",
-    # assuming mu=5e-10 (from Drake 1991).
-    population_size=1.8e8,
-    population_size_citations=[_hartl_et_al.because(stdpopsim.CiteReason.POP_SIZE)],
-)
-
-
-stdpopsim.register_species(_species)
+from . import species  # noqa 401
diff --git a/stdpopsim/catalog/EscCol/species.py b/stdpopsim/catalog/EscCol/species.py
new file mode 100644
index 000000000..a95082c8c
--- /dev/null
+++ b/stdpopsim/catalog/EscCol/species.py
@@ -0,0 +1,76 @@
+import stdpopsim
+from . import genome_data
+
+###########################################################
+#
+# Genome definition
+#
+###########################################################
+
+_hartl_et_al = stdpopsim.Citation(
+    author="Hartl, Moriyama, and Sawyer",
+    year="1994",
+    # doesn't have a doi
+    doi="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1206133/",
+)
+
+_sezonov_et_al = stdpopsim.Citation(
+    author="Sezonov et al.", year="2007", doi="https://doi.org/10.1128/JB.01368-07"
+)
+
+_wielgoss_et_al = stdpopsim.Citation(
+    author="Wielgoss et al.", year="2011", doi="https://doi.org/10.1534/g3.111.000406"
+)
+
+_blattner_et_al = stdpopsim.Citation(
+    author="Blattner et al.", year="1997", doi="10.1126/science.277.5331.1453"
+)
+
+_chromosomes = []
+for name, data in genome_data.data["chromosomes"].items():
+    _chromosomes.append(
+        stdpopsim.Chromosome(
+            id=name,
+            length=data["length"],
+            synonyms=data["synonyms"],
+            # Wielgoss et al. (2011) calculated for strain REL606,
+            # from synonymous substitutions over 40,000 generations.
+            mutation_rate=8.9e-11,
+            recombination_rate=0.0,
+        )
+    )
+
+# mean_conversion_rate=8.9e-11 # not implemented yet!
+# mean_conversion_length=542 # not implemented yet!
+
+_genome = stdpopsim.Genome(
+    chromosomes=_chromosomes,
+    assembly_name=genome_data.data["assembly_name"],
+    assembly_accession=genome_data.data["assembly_accession"],
+    mutation_rate_citations=[
+        _wielgoss_et_al.because(stdpopsim.CiteReason.MUT_RATE),
+    ],
+    assembly_citations=[_blattner_et_al.because(stdpopsim.CiteReason.ASSEMBLY)],
+)
+
+
+_species = stdpopsim.Species(
+    id="EscCol",
+    name="Escherichia coli",
+    common_name="E. coli",
+    # We use the K-12 strain, because the parameters we're using more
+    # closely match this strain than the ensembl default (HUSEC2011).
+    ensembl_id="escherichia_coli_str_k_12_substr_mg1655",
+    genome=_genome,
+    # E. coli K-12 strain MG1655 "doubling time during steady-state growth in
+    # Luria-Bertani broth was 20 min".
+    generation_time=0.00003805175,  # 1.0 / (525600 min/year / 20 min/gen)
+    generation_time_citations=[_sezonov_et_al.because(stdpopsim.CiteReason.GEN_TIME)],
+    # Hartl et al. calculated Ne for "natural isolates of E. coli",
+    # assuming mu=5e-10 (from Drake 1991).
+    population_size=1.8e8,
+    population_size_citations=[_hartl_et_al.because(stdpopsim.CiteReason.POP_SIZE)],
+)
+
+
+stdpopsim.register_species(_species)
diff --git a/stdpopsim/catalog/HomSap/__init__.py b/stdpopsim/catalog/HomSap/__init__.py
index 18499ec7a..f20451bd9 100644
--- a/stdpopsim/catalog/HomSap/__init__.py
+++ b/stdpopsim/catalog/HomSap/__init__.py
@@ -1,1807 +1,7 @@
 """
 Catalog definitions for Homo Sapiens
 """
-import math
-
-import msprime
-
-import stdpopsim
-from . import genome_data
-
-###########################################################
-#
-# Genome definition
-#
-###########################################################
-
-
-# mean_recombination_rate was computed across all windows of the HapMapII_GRCh37 map
-_recombination_rate_data = {
-    "1": 1.1485597641285933e-08,
-    "2": 1.1054289277533446e-08,
-    "3": 1.1279585624662551e-08,
-    "4": 1.1231162636001008e-08,
-    "5": 1.1280936570022824e-08,
-    "6": 1.1222852661225285e-08,
-    "7": 1.1764614397655721e-08,
-    "8": 1.1478465778920576e-08,
-    "9": 1.1780701596308656e-08,
-    "10": 1.3365134257075317e-08,
-    "11": 1.1719334320833283e-08,
-    "12": 1.305017186986983e-08,
-    "13": 1.0914860554958317e-08,
-    "14": 1.119730771394731e-08,
-    "15": 1.3835785893339787e-08,
-    "16": 1.4834607113882717e-08,
-    "17": 1.582489036239487e-08,
-    "18": 1.5075956950023575e-08,
-    "19": 1.8220141872466202e-08,
-    "20": 1.7178269031631664e-08,
-    "21": 1.3045214034879191e-08,
-    "22": 1.4445022767788226e-08,
-    "X": 1.164662223273842e-08,
-    "Y": 0.0,
-    "MT": 0.0,
-}
-
-_genome2001 = stdpopsim.Citation(
-    doi="http://dx.doi.org/10.1038/35057062",
-    year="2001",
-    author="International Human Genome Sequencing Consortium",
-    reasons={stdpopsim.CiteReason.ASSEMBLY},
-)
-
-_hapmap2007 = stdpopsim.Citation(
-    doi="https://doi.org/10.1038/nature06258",
-    year=2007,
-    author="The International HapMap Consortium",
-)
-
-_takahata1993 = stdpopsim.Citation(
-    doi="https://doi.org/10.1093/oxfordjournals.molbev.a039995",
-    year="1993",
-    author="Takahata",
-    reasons={stdpopsim.CiteReason.POP_SIZE},
-)
-
-_tian2019 = stdpopsim.Citation(
-    doi="https://doi.org/10.1016/j.ajhg.2019.09.012",
-    year="2019",
-    author="Tian, Browning, and Browning",
-    reasons={stdpopsim.CiteReason.MUT_RATE},
-)
-
-_tremblay2000 = stdpopsim.Citation(
-    doi="https://doi.org/10.1086/302770",
-    year="2000",
-    author="Tremblay and Vézina",
-    reasons={stdpopsim.CiteReason.GEN_TIME},
-)
-
-_chromosomes = []
-for name, data in genome_data.data["chromosomes"].items():
-    _chromosomes.append(
-        stdpopsim.Chromosome(
-            id=name,
-            length=data["length"],
-            synonyms=data["synonyms"],
-            mutation_rate=1.29e-8,
-            recombination_rate=_recombination_rate_data[name],
-        )
-    )
-
-_genome = stdpopsim.Genome(
-    chromosomes=_chromosomes,
-    mutation_rate_citations=[_tian2019.because(stdpopsim.CiteReason.MUT_RATE)],
-    recombination_rate_citations=[_hapmap2007.because(stdpopsim.CiteReason.REC_RATE)],
-    assembly_name=genome_data.data["assembly_name"],
-    assembly_accession=genome_data.data["assembly_accession"],
-    assembly_citations=[_genome2001],
-)
-
-_species = stdpopsim.Species(
-    id="HomSap",
-    name="Homo sapiens",
-    common_name="Human",
-    genome=_genome,
-    generation_time=30,
-    generation_time_citations=[_tremblay2000.because(stdpopsim.CiteReason.GEN_TIME)],
-    population_size=10 ** 4,
-    population_size_citations=[_takahata1993.because(stdpopsim.CiteReason.POP_SIZE)],
-)
-
-stdpopsim.register_species(_species)
-
-
-###########################################################
-#
-# Genetic maps
-#
-###########################################################
-
-
-_gm = stdpopsim.GeneticMap(
-    species=_species,
-    id="HapMapII_GRCh37",
-    description="HapMap Phase II lifted over to GRCh37",
-    long_description="""
-        This genetic map is from the Phase II Hapmap project
-        and based on 3.1 million genotyped SNPs
-        from 270 individuals across four populations (YRI, CEU, CHB and JPT).
-        Genome wide recombination rates were estimated using LDHat.
-        This version of the HapMap genetic map was lifted over to GRCh37
-        (and adjusted in regions where the genome assembly had rearranged)
-        for use in the 1000 Genomes project. Please see the README file on
-        the 1000 Genomes download site for details of these adjustments.
-        ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/working/20110106_recombination_hotspots
-        """,
-    url=(
-        "https://stdpopsim.s3-us-west-2.amazonaws.com/genetic_maps/"
-        "HomSap/HapmapII_GRCh37_RecombinationHotspots.tar.gz"
-    ),
-    sha256="80f22d9e6cb0e497074ed1bc277e765fa9d8e22f21b2f66c3b10286520f6b68f",
-    file_pattern="genetic_map_GRCh37_chr{id}.txt",
-    citations=[_hapmap2007.because(stdpopsim.CiteReason.GEN_MAP)],
-)
-_species.add_genetic_map(_gm)
-
-_gm = stdpopsim.GeneticMap(
-    species=_species,
-    id="DeCodeSexAveraged_GRCh36",
-    description="Sex averaged map from deCode family study",
-    long_description="""
-        This genetic map is from the deCode study of recombination
-        events in 15,257 parent-offspring pairs from Iceland.
-        289,658 phased autosomal SNPs were used to call recombinations
-        within these families, and recombination rates computed from the
-        density of these events. This is the combined male and female
-        (sex averaged) map. See
-        https://www.decode.com/addendum/ for more details.""",
-    url=(
-        "https://stdpopsim.s3-us-west-2.amazonaws.com/genetic_maps/"
-        "HomSap/decode_2010_sex-averaged_map.tar.gz"
-    ),
-    sha256="a17f491756f863c1cc61240c0ce4f383d8085c68b45edd3fed677bb9692af23d",
-    file_pattern="genetic_map_decode_2010_sex-averaged_chr{id}.txt",
-    citations=[
-        stdpopsim.Citation(
-            year=2010,
-            author="Kong et al",
-            doi="https://doi.org/10.1038/nature09525",
-            reasons={stdpopsim.CiteReason.GEN_MAP},
-        )
-    ],
-)
-_species.add_genetic_map(_gm)
-
-# 26 populations in 1000 Genomes
-for pop, sha256 in zip(
-    [
-        "ACB",
-        "ASW",
-        "BEB",
-        "CDX",
-        "CEU",
-        "CHB",
-        "CHS",
-        "CLM",
-        "ESN",
-        "FIN",
-        "GBR",
-        "GIH",
-        "GWD",
-        "IBS",
-        "ITU",
-        "JPT",
-        "KHV",
-        "LWK",
-        "MSL",
-        "MXL",
-        "PEL",
-        "PJL",
-        "PUR",
-        "STU",
-        "TSI",
-        "YRI",
-    ],
-    """
-588be112b6aac9f82bfb64ea828809051cdf021a627450b2d8c37d7909486903
-e7017a82e92e87f443cd21c7f1c3ce24168fa5efaa48f8448c8e095d13aebcaf
-c6f69e6a847af6f88fbf1c82024af6fd5eeb35c167e3a721d9de6d6ff146a3ba
-37cd4e8ac16775b9c8126dc09ab90840925f47138675d6a5942dd6264d089f6c
-fcfacc1858bff5a64480982bab0821d47c4c56262310dd85290f907b8ea6b081
-3836a7853b2e128b6bdbeb0e5eb1ba582975304a9efa10b742be856efdc520a9
-e0ba174edf944e96f327d799010a5be8ecc110968be6b27fb8f6901d5a413cbb
-969f04541d90b70e68f75c0edcfffeb00bee2cd97cf7ad598439f18a8d3a0885
-6310f3860e4ccfcedddc173df2a03400d654ddd83c78eb88ec92e104d00bf40b
-bb7ac5df072b8b2a4567825c5ee4382b67cae5290f5879a1d7c3f743f244b3f6
-798daacb41ae538017303f103c39cdb88add7d88eb964922fe6c81406d8aaa64
-28e4552e485b55b0bc392485cd389a155fa83ee3556f221b3b3377d2c996cbdf
-f2d57c835ac61763112c15acd9cc8240af7c7fb1d7cfa7543cb65a565da2f62b
-bbc49af3646b47dd5b58b8ae28716dcce695eb4a9cb223e261f45345b1a051b9
-050b965efb0b4c9197a289521406586f4207dd4833be9fd429d118f7d061788a
-e7284218b4bf3e1487b9194c9f243c5dff88f5edbd40247118225a9acc08c76e
-74cd0bc6bc12c2b59f928c6a3d60773b6118625375266ffba0c89da61c376f55
-dd3ef07641ad7c8d3bfa74a98a5820f1b70ea74c3e82103663bbe61931bf8579
-ad5064ee575a10f890c81b9bdbd55128197582db1c9b153b575f5e7fbe542b0d
-f37b50b940816d9c12429a4b71709b18569e77c1ea4c78d5fed19141751822c3
-024ac242da742d006a210cbf77c5b8457735a6690dc78cc01686deafe7cb7beb
-871b518e3b9ca92838e206125e4a1187690018ec5bddeaf39eb2e6f6be5705f1
-00f87ca1ad6430b272251d3ce5fa93b44a0a50d6472237bcb5dbad2dbcd567b8
-fa9e61492de4c160cf987ac6a8cc3a77efdea751f0f5aa9303a294a37f81cbd7
-e3cf06041c6ffceb29fbc617c86966948203e166c3399d036c95915e3d5ebdcd
-2989aecaa086af37cdef9ee2245f12000b3f29c08b5d1a71744578a684fe377f
-    """.strip().splitlines(),
-):
-    _gm = stdpopsim.GeneticMap(
-        species=_species,
-        id="Pyrho{}_GRCh38".format(pop),
-        description="Pyrho population-specific map for {}".format(pop),
-        long_description="""
-        This genetic map was inferred using individuals from the {}
-        population from Phase 3 of the 1000 Genomes Project. Rates
-        were estimated using pyrho (https://github.com/popgenmethods/pyrho)
-        while using population-specific population size history estimates
-        obtained from smc++ (https://github.com/popgenmethods/smcpp).
-        Genetic maps are only available for the 22 autosomes.
-        See https://doi.org/10.1126/sciadv.aaw9206 for more
-        details.""".format(
-            pop
-        ),
-        url=(
-            "https://stdpopsim.s3-us-west-2.amazonaws.com/genetic_maps/"
-            "HomSap/Pyrho{}_GRCh38.tar.gz".format(pop)
-        ),
-        sha256=sha256,
-        file_pattern=("Pyrho{}_GRCh38_chr{{id}}.txt".format(pop)),
-        citations=[
-            stdpopsim.Citation(
-                year=2019,
-                author="Spence and Song",
-                doi="https://doi.org/10.1126/sciadv.aaw9206",
-                reasons={stdpopsim.CiteReason.GEN_MAP},
-            )
-        ],
-    )
-    _species.add_genetic_map(_gm)
-
-###########################################################
-#
-# Annotations
-#
-#
-#
-###########################################################
-
-_an = stdpopsim.Annotation(
-    species=_species,
-    id="Ensembl_GRCh38_gff3",
-    description="Ensembl GFF3 annotations on GRCh38",
-    url=(
-        "ftp://ftp.ensembl.org/pub/release-101/"
-        "gff3/homo_sapiens/Homo_sapiens.GRCh38.101.gff3.gz"
-    ),
-    zarr_url=(
-        "https://stdpopsim.s3-us-west-2.amazonaws.com/annotations/HomSap.GRCh38.zip"
-    ),
-    zarr_sha256="929c52dc5e5172d4a96a776d066f014f8f207f1477a21f7b88626f81a7142d0b",
-    citations=[
-        stdpopsim.Citation(
-            year=2018,
-            author="Hunt et al",
-            doi="https://doi.org/10.1093/database/bay119",
-            reasons={stdpopsim.CiteReason.ANNOTATION},
-        )
-    ],
-)
-_species.add_annotations(_an)
-
-###########################################################
-#
-# Demographic models
-#
-###########################################################
-
-# population definitions that are reused.
-_yri_population = stdpopsim.Population(
-    id="YRI", description="1000 Genomes YRI (Yorubans)"
-)
-_ceu_population = stdpopsim.Population(
-    id="CEU",
-    description=(
-        "1000 Genomes CEU (Utah Residents (CEPH) with Northern and "
-        "Western European Ancestry"
-    ),
-)
-_chb_population = stdpopsim.Population(
-    id="CHB", description="1000 Genomes CHB (Han Chinese in Beijing, China)"
-)
-
-
-_tennessen_et_al = stdpopsim.Citation(
-    author="Tennessen et al.",
-    year=2012,
-    doi="https://doi.org/10.1126/science.1219240",
-    reasons={stdpopsim.CiteReason.DEM_MODEL},
-)
-
-
-def _ooa_3():
-    id = "OutOfAfrica_3G09"
-    description = "Three population out-of-Africa"
-    long_description = """
-        The three population Out-of-Africa model from Gutenkunst et al. 2009.
-        It describes the ancestral human population in Africa, the out of Africa
-        event, and the subsequent European-Asian population split.
-        Model parameters are the maximum likelihood values of the
-        various parameters given in Table 1 of Gutenkunst et al.
-    """
-    populations = [_yri_population, _ceu_population, _chb_population]
-
-    citations = [
-        stdpopsim.Citation(
-            author="Gutenkunst et al.",
-            year=2009,
-            doi="https://doi.org/10.1371/journal.pgen.1000695",
-            reasons={stdpopsim.CiteReason.DEM_MODEL},
-        )
-    ]
-
-    generation_time = 25
-
-    # First we set out the maximum likelihood values of the various parameters
-    # given in Table 1.
-    N_A = 7300
-    N_B = 2100
-    N_AF = 12300
-    N_EU0 = 1000
-    N_AS0 = 510
-    # Times are provided in years, so we convert into generations.
-
-    T_AF = 220e3 / generation_time
-    T_B = 140e3 / generation_time
-    T_EU_AS = 21.2e3 / generation_time
-    # We need to work out the starting (diploid) population sizes based on
-    # the growth rates provided for these two populations
-    r_EU = 0.004
-    r_AS = 0.0055
-    N_EU = N_EU0 / math.exp(-r_EU * T_EU_AS)
-    N_AS = N_AS0 / math.exp(-r_AS * T_EU_AS)
-    # Migration rates during the various epochs.
-    m_AF_B = 25e-5
-    m_AF_EU = 3e-5
-    m_AF_AS = 1.9e-5
-    m_EU_AS = 9.6e-5
-
-    return stdpopsim.DemographicModel(
-        id=id,
-        description=description,
-        long_description=long_description,
-        populations=populations,
-        citations=citations,
-        generation_time=generation_time,
-        # Population IDs correspond to their indexes in the population
-        # configuration array. Therefore, we have 0=YRI, 1=CEU and 2=CHB
-        # initially.
-        population_configurations=[
-            msprime.PopulationConfiguration(
-                initial_size=N_AF, metadata=populations[0].asdict()
-            ),
-            msprime.PopulationConfiguration(
-                initial_size=N_EU, growth_rate=r_EU, metadata=populations[1].asdict()
-            ),
-            msprime.PopulationConfiguration(
-                initial_size=N_AS, growth_rate=r_AS, metadata=populations[2].asdict()
-            ),
-        ],
-        migration_matrix=[
-            [0, m_AF_EU, m_AF_AS],  # noqa
-            [m_AF_EU, 0, m_EU_AS],  # noqa
-            [m_AF_AS, m_EU_AS, 0],  # noqa
-        ],
-        demographic_events=[
-            # CEU and CHB merge into B with rate changes at T_EU_AS
-            msprime.MassMigration(
-                time=T_EU_AS, source=2, destination=1, proportion=1.0
-            ),
-            msprime.MigrationRateChange(time=T_EU_AS, rate=0),
-            msprime.MigrationRateChange(time=T_EU_AS, rate=m_AF_B, matrix_index=(0, 1)),
-            msprime.MigrationRateChange(time=T_EU_AS, rate=m_AF_B, matrix_index=(1, 0)),
-            msprime.PopulationParametersChange(
-                time=T_EU_AS, initial_size=N_B, growth_rate=0, population_id=1
-            ),
-            # Population B merges into YRI at T_B
-            msprime.MassMigration(time=T_B, source=1, destination=0, proportion=1.0),
-            msprime.MigrationRateChange(time=T_B, rate=0),
-            # Size changes to N_A at T_AF
-            msprime.PopulationParametersChange(
-                time=T_AF, initial_size=N_A, population_id=0
-            ),
-        ],
-    )
-
-
-_species.add_demographic_model(_ooa_3())
-
-
-def _ooa_2():
-    id = "OutOfAfrica_2T12"
-    description = "Two population out-of-Africa"
-    long_description = """
-        The model is derived from the Tennesen et al. analysis of the
-        jSFS from European Americans and African Americans.
-        It describes the ancestral human population in Africa, the out of Africa event,
-        and two distinct periods of subsequent European population growth over the past
-        23kya. Model parameters are taken from Fig. S5 in Fu et al.
-    """
-    populations = [
-        stdpopsim.Population(id="AFR", description="African Americans"),
-        stdpopsim.Population(id="EUR", description="European Americans"),
-    ]
-    citations = [
-        _tennessen_et_al,
-        stdpopsim.Citation(
-            author="Fu et al.",
-            year=2013,
-            doi="https://doi.org/10.1038/nature11690",
-            reasons={stdpopsim.CiteReason.DEM_MODEL},
-        ),
-    ]
-
-    generation_time = 25
-
-    T_AF = 148e3 / generation_time
-    T_OOA = 51e3 / generation_time
-    T_EU0 = 23e3 / generation_time
-    T_EG = 5115 / generation_time
-
-    # Growth rates
-    r_EU0 = 0.00307
-    r_EU = 0.0195
-    r_AF = 0.0166
-
-    # population sizes
-    N_A = 7310
-    N_AF1 = 14474
-    N_B = 1861
-    N_EU0 = 1032
-    N_EU1 = N_EU0 / math.exp(-r_EU0 * (T_EU0 - T_EG))
-
-    # migration rates
-    m_AF_B = 15e-5
-    m_AF_EU = 2.5e-5
-
-    # present Ne
-    N_EU = N_EU1 / math.exp(-r_EU * T_EG)
-    N_AF = N_AF1 / math.exp(-r_AF * T_EG)
-
-    return stdpopsim.DemographicModel(
-        id=id,
-        description=description,
-        long_description=long_description,
-        populations=populations,
-        citations=citations,
-        generation_time=generation_time,
-        population_configurations=[
-            msprime.PopulationConfiguration(
-                initial_size=N_AF, growth_rate=r_AF, metadata=populations[0].asdict()
-            ),
-            msprime.PopulationConfiguration(
-                initial_size=N_EU, growth_rate=r_EU, metadata=populations[1].asdict()
-            ),
-        ],
-        migration_matrix=[
-            [0, m_AF_EU],
-            [m_AF_EU, 0],
-        ],
-        demographic_events=[
-            msprime.MigrationRateChange(time=T_EG, rate=m_AF_EU, matrix_index=(0, 1)),
-            msprime.MigrationRateChange(time=T_EG, rate=m_AF_EU, matrix_index=(1, 0)),
-            msprime.PopulationParametersChange(
-                time=T_EG, growth_rate=r_EU0, initial_size=N_EU1, population_id=1
-            ),
-            msprime.PopulationParametersChange(
-                time=T_EG, growth_rate=0, initial_size=N_AF1, population_id=0
-            ),
-            msprime.MigrationRateChange(time=T_EU0, rate=m_AF_B, matrix_index=(0, 1)),
-            msprime.MigrationRateChange(time=T_EU0, rate=m_AF_B, matrix_index=(1, 0)),
-            msprime.PopulationParametersChange(
-                time=T_EU0, initial_size=N_B, growth_rate=0, population_id=1
-            ),
-            msprime.MassMigration(time=T_OOA, source=1, destination=0, proportion=1.0),
-            msprime.MigrationRateChange(time=T_OOA, rate=0),
-            msprime.PopulationParametersChange(
-                time=T_AF, initial_size=N_A, population_id=0
-            ),
-        ],
-    )
-
-
-_species.add_demographic_model(_ooa_2())
-
-
-def _african():
-    id = "Africa_1T12"
-    description = "African population"
-    long_description = """
-        The model is a simplification of the two population Tennesen et al.
-        model with the European-American population removed so that we are
-        modeling the African population in isolation.
-    """
-    populations = [
-        stdpopsim.Population(id="AFR", description="African"),
-    ]
-    citations = [_tennessen_et_al]
-
-    generation_time = 25
-
-    T_AF = 148e3 / generation_time
-    T_EG = 5115 / generation_time
-
-    # Growth rate
-    r_AF = 0.0166
-
-    # population sizes
-    N_A = 7310
-    N_AF1 = 14474
-
-    # present Ne
-    N_AF = N_AF1 / math.exp(-r_AF * T_EG)
-
-    return stdpopsim.DemographicModel(
-        id=id,
-        description=description,
-        long_description=long_description,
-        populations=populations,
-        citations=citations,
-        generation_time=generation_time,
-        population_configurations=[
-            msprime.PopulationConfiguration(
-                initial_size=N_AF, growth_rate=r_AF, metadata=populations[0].asdict()
-            ),
-        ],
-        demographic_events=[
-            msprime.PopulationParametersChange(
-                time=T_EG, growth_rate=0, initial_size=N_AF1, population_id=0
-            ),
-            msprime.PopulationParametersChange(
-                time=T_AF, initial_size=N_A, population_id=0
-            ),
-        ],
-    )
-
-
-_species.add_demographic_model(_african())
-
-
-def _america():
-    id = "AmericanAdmixture_4B11"
-    description = "American admixture"
-    long_description = """
-        Demographic model for American admixture, taken from Browning et al. 2011.
-        This model extends the Gravel et al. (2011) model of African/European/Asian
-        demographic history to simulate an admixed population with admixture
-        occurring 12 generations ago. The admixed population had an initial size
-        of 30,000 and grew at a rate of 5% per generation, with 1/6 of the
-        population of African ancestry, 1/3 European, and 1/2 Asian.
-    """
-    populations = [
-        stdpopsim.Population(id="AFR", description="Contemporary African population"),
-        stdpopsim.Population(id="EUR", description="Contemporary European population"),
-        stdpopsim.Population(id="ASIA", description="Contemporary Asian population"),
-        stdpopsim.Population(id="ADMIX", description="Modern admixed population"),
-    ]
-
-    citations = [
-        stdpopsim.Citation(
-            author="Browning et al.",
-            year=2011,
-            doi="http://dx.doi.org/10.1371/journal.pgen.1007385",
-            reasons={stdpopsim.CiteReason.DEM_MODEL},
-        )
-    ]
-
-    generation_time = 25
-
-    # Model code was ported from Supplementary File 1.
-    N0 = 7310  # initial population size
-    Thum = 5920  # time (gens) of advent of modern humans
-    Naf = 14474  # size of african population
-    Tooa = 2040  # number of generations back to Out of Africa
-    Nb = 1861  # size of out of Africa population
-    mafb = 1.5e-4  # migration rate Africa and Out-of-Africa
-    Teu = 920  # number generations back to Asia-Europe split
-    Neu = 1032  # bottleneck population sizes
-    Nas = 554
-    mafeu = 2.5e-5  # mig. rates
-    mafas = 7.8e-6
-    meuas = 3.11e-5
-    reu = 0.0038  # growth rate per generation in Europe
-    ras = 0.0048  # growth rate per generation in Asia
-    Tadmix = 12  # time of admixture
-    Nadmix = 30000  # initial size of admixed population
-    radmix = 0.05  # growth rate of admixed population
-    # pop0 is Africa, pop1 is Europe, pop2 is Asia,  pop3 is admixed
-
-    population_configurations = [
-        msprime.PopulationConfiguration(
-            initial_size=Naf, growth_rate=0.0, metadata=populations[0].asdict()
-        ),
-        msprime.PopulationConfiguration(
-            initial_size=Neu * math.exp(reu * Teu),
-            growth_rate=reu,
-            metadata=populations[1].asdict(),
-        ),
-        msprime.PopulationConfiguration(
-            initial_size=Nas * math.exp(ras * Teu),
-            growth_rate=ras,
-            metadata=populations[2].asdict(),
-        ),
-        msprime.PopulationConfiguration(
-            initial_size=Nadmix * math.exp(radmix * Tadmix),
-            growth_rate=radmix,
-            metadata=populations[3].asdict(),
-        ),
-    ]
-
-    migration_matrix = [
-        [0, mafeu, mafas, 0],
-        [mafeu, 0, meuas, 0],
-        [mafas, meuas, 0, 0],
-        [0, 0, 0, 0],
-    ]
-    # Admixture event, 1/6 Africa, 2/6 Europe, 3/6 Asia
-    admixture_event = [
-        msprime.MassMigration(
-            time=Tadmix, source=3, destination=0, proportion=1.0 / 6.0
-        ),
-        msprime.MassMigration(
-            time=Tadmix, source=3, destination=1, proportion=2.0 / 5.0
-        ),
-        msprime.MassMigration(time=Tadmix, source=3, destination=2, proportion=1.0),
-    ]
-    # Asia and Europe split
-    eu_event = [
-        msprime.MigrationRateChange(time=Teu, rate=0.0),
-        msprime.MassMigration(time=Teu, source=2, destination=1, proportion=1.0),
-        msprime.PopulationParametersChange(
-            time=Teu, initial_size=Nb, growth_rate=0.0, population_id=1
-        ),
-        msprime.MigrationRateChange(time=Teu, rate=mafb, matrix_index=(0, 1)),
-        msprime.MigrationRateChange(time=Teu, rate=mafb, matrix_index=(1, 0)),
-    ]
-    # Out of Africa event
-    ooa_event = [
-        msprime.MigrationRateChange(time=Tooa, rate=0.0),
-        msprime.MassMigration(time=Tooa, source=1, destination=0, proportion=1.0),
-    ]
-    # initial population size
-    init_event = [
-        msprime.PopulationParametersChange(time=Thum, initial_size=N0, population_id=0)
-    ]
-    demographic_events = admixture_event + eu_event + ooa_event + init_event
-
-    return stdpopsim.DemographicModel(
-        id=id,
-        description=description,
-        long_description=long_description,
-        populations=populations,
-        citations=citations,
-        generation_time=generation_time,
-        population_configurations=population_configurations,
-        migration_matrix=migration_matrix,
-        demographic_events=demographic_events,
-    )
-
-
-_species.add_demographic_model(_america())
-
-
-def _ooa_archaic():
-    id = "OutOfAfricaArchaicAdmixture_5R19"
-    description = "Three population out-of-Africa with archaic admixture"
-    long_description = """
-        The three population out-of-African model popularized by Gutenkunst et al. (2009)
-        and augmented by archaic contributions to both Eurasian and African populations.
-        Two archaic populations split early in human history, before the African
-        expansion, and contribute to Eurasian populations (putative Neanderthal branch)
-        and to the African branch (a deep diverging branch within Africa). Admixture
-        is modeled as symmetric migration between the archaic and modern human branches,
-        with contribution ending at a given time in the past.
-    """
-    populations = [
-        _yri_population,
-        _ceu_population,
-        _chb_population,
-        stdpopsim.Population(
-            "Neanderthal", "Putative Neanderthals", sampling_time=None
-        ),
-        stdpopsim.Population(
-            "ArchaicAFR", "Putative Archaic Africans", sampling_time=None
-        ),
-    ]
-    citations = [
-        stdpopsim.Citation(
-            author="Ragsdale and Gravel",
-            year=2019,
-            doi="https://doi.org/10.1371/journal.pgen.1008204",
-            reasons={stdpopsim.CiteReason.DEM_MODEL},
-        )
-    ]
-
-    # First we set out the maximum likelihood values of the various parameters
-    # given in Table 1 (under archaic admixture).
-    N_0 = 3600
-    N_YRI = 13900
-    N_B = 880
-    N_CEU0 = 2300
-    N_CHB0 = 650
-
-    # Times are provided in years, so we convert into generations.
-    # In the published model, the authors used a generation time of 29 years to
-    # convert from genetic to physical units
-    generation_time = 29
-
-    T_AF = 300e3 / generation_time
-    T_B = 60.7e3 / generation_time
-    T_EU_AS = 36.0e3 / generation_time
-    T_arch_afr_split = 499e3 / generation_time
-    T_arch_afr_mig = 125e3 / generation_time
-    T_nean_split = 559e3 / generation_time
-    T_arch_adm_end = 18.7e3 / generation_time
-
-    # We need to work out the starting (diploid) population sizes based on
-    # the growth rates provided for these two populations
-    r_CEU = 0.00125
-    r_CHB = 0.00372
-    N_CEU = N_CEU0 / math.exp(-r_CEU * T_EU_AS)
-    N_CHB = N_CHB0 / math.exp(-r_CHB * T_EU_AS)
-
-    # Migration rates during the various epochs.
-    m_AF_B = 52.2e-5
-    m_YRI_CEU = 2.48e-5
-    m_YRI_CHB = 0e-5
-    m_CEU_CHB = 11.3e-5
-    m_AF_arch_af = 1.98e-5
-    m_OOA_nean = 0.825e-5
-
-    # Population IDs correspond to their indexes in the population
-    # configuration array. Therefore, we have 0=YRI, 1=CEU and 2=CHB
-    # initially.
-    # We also have two archaic populations, putative Neanderthals and
-    # archaicAfrican, which are population indices 3=Nean and 4=arch_afr.
-    # Their sizes are equal to the ancestral reference population size N_0.
-    population_configurations = [
-        msprime.PopulationConfiguration(
-            initial_size=N_YRI, metadata=populations[0].asdict()
-        ),
-        msprime.PopulationConfiguration(
-            initial_size=N_CEU, growth_rate=r_CEU, metadata=populations[1].asdict()
-        ),
-        msprime.PopulationConfiguration(
-            initial_size=N_CHB, growth_rate=r_CHB, metadata=populations[2].asdict()
-        ),
-        msprime.PopulationConfiguration(
-            initial_size=N_0, metadata=populations[3].asdict()
-        ),
-        msprime.PopulationConfiguration(
-            initial_size=N_0, metadata=populations[4].asdict()
-        ),
-    ]
-    migration_matrix = [  # noqa
-        [0, m_YRI_CEU, m_YRI_CHB, 0, 0],  # noqa
-        [m_YRI_CEU, 0, m_CEU_CHB, 0, 0],  # noqa
-        [m_YRI_CHB, m_CEU_CHB, 0, 0, 0],  # noqa
-        [0, 0, 0, 0, 0],  # noqa
-        [0, 0, 0, 0, 0],  # noqa
-    ]  # noqa
-    demographic_events = [
-        # first event is migration turned on between modern and archaic humans
-        msprime.MigrationRateChange(
-            time=T_arch_adm_end, rate=m_AF_arch_af, matrix_index=(0, 4)
-        ),
-        msprime.MigrationRateChange(
-            time=T_arch_adm_end, rate=m_AF_arch_af, matrix_index=(4, 0)
-        ),
-        msprime.MigrationRateChange(
-            time=T_arch_adm_end, rate=m_OOA_nean, matrix_index=(1, 3)
-        ),
-        msprime.MigrationRateChange(
-            time=T_arch_adm_end, rate=m_OOA_nean, matrix_index=(3, 1)
-        ),
-        msprime.MigrationRateChange(
-            time=T_arch_adm_end, rate=m_OOA_nean, matrix_index=(2, 3)
-        ),
-        msprime.MigrationRateChange(
-            time=T_arch_adm_end, rate=m_OOA_nean, matrix_index=(3, 2)
-        ),
-        # CEU and CHB merge into B with rate changes at T_EU_AS
-        msprime.MassMigration(time=T_EU_AS, source=2, destination=1, proportion=1.0),
-        msprime.MigrationRateChange(time=T_EU_AS, rate=0),
-        msprime.MigrationRateChange(time=T_EU_AS, rate=m_AF_B, matrix_index=(0, 1)),
-        msprime.MigrationRateChange(time=T_EU_AS, rate=m_AF_B, matrix_index=(1, 0)),
-        msprime.MigrationRateChange(
-            time=T_EU_AS, rate=m_AF_arch_af, matrix_index=(0, 4)
-        ),
-        msprime.MigrationRateChange(
-            time=T_EU_AS, rate=m_AF_arch_af, matrix_index=(4, 0)
-        ),
-        msprime.MigrationRateChange(time=T_EU_AS, rate=m_OOA_nean, matrix_index=(1, 3)),
-        msprime.MigrationRateChange(time=T_EU_AS, rate=m_OOA_nean, matrix_index=(3, 1)),
-        msprime.PopulationParametersChange(
-            time=T_EU_AS, initial_size=N_B, growth_rate=0, population_id=1
-        ),
-        # Population B merges into YRI at T_B
-        msprime.MassMigration(time=T_B, source=1, destination=0, proportion=1.0),
-        msprime.MigrationRateChange(time=T_B, rate=0),
-        msprime.MigrationRateChange(time=T_B, rate=m_AF_arch_af, matrix_index=(0, 4)),
-        msprime.MigrationRateChange(time=T_B, rate=m_AF_arch_af, matrix_index=(4, 0)),
-        # Beginning of migration between African and archaic African populations
-        msprime.MigrationRateChange(time=T_arch_afr_mig, rate=0),
-        # Size changes to N_0 at T_AF
-        msprime.PopulationParametersChange(
-            time=T_AF, initial_size=N_0, population_id=0
-        ),
-        # Archaic African merges with moderns
-        msprime.MassMigration(
-            time=T_arch_afr_split, source=4, destination=0, proportion=1.0
-        ),
-        # Neanderthal merges with moderns
-        msprime.MassMigration(
-            time=T_nean_split, source=3, destination=0, proportion=1.0
-        ),
-    ]
-
-    return stdpopsim.DemographicModel(
-        id=id,
-        description=description,
-        long_description=long_description,
-        populations=populations,
-        citations=citations,
-        generation_time=generation_time,
-        population_configurations=population_configurations,
-        migration_matrix=migration_matrix,
-        demographic_events=demographic_events,
-    )
-
-
-_species.add_demographic_model(_ooa_archaic())
-
-
-def _zigzag():
-    id = "Zigzag_1S14"
-    description = "Periodic growth and decline."
-    long_description = """
-        A validation model used by Schiffels and Durbin (2014) and Terhorst and
-        Terhorst, Kamm, and Song (2017) with periods of exponential growth and
-        decline in a single population.
-        """
-    populations = [
-        stdpopsim.Population("generic", "Generic expanding and contracting population"),
-    ]
-    citations = [
-        stdpopsim.Citation(
-            author="Schiffels and Durbin",
-            year=2014,
-            doi="https://doi.org/10.1038/ng.3015",
-            reasons={stdpopsim.CiteReason.DEM_MODEL},
-        )
-    ]
-
-    generation_time = 30
-    N0 = 5 * 14312
-
-    g_1 = 0.02302578256
-    t_1 = 33.333334976
-
-    g_2 = -0.0057564631
-    t_2 = 133.3334544
-
-    g_3 = 0.00143911578
-    t_3 = 533.33324512
-
-    g_4 = -0.0003597785
-    t_4 = 2133.3352704
-
-    g_5 = 8.9944801565e-5
-    t_5 = 8533.329632
-
-    n_ancient = N0 / 10
-    t_ancient = 34133.318528
-
-    population_configurations = [
-        msprime.PopulationConfiguration(
-            initial_size=N0, metadata=populations[0].asdict()
-        )
-    ]
-
-    demographic_events = [
-        msprime.PopulationParametersChange(time=t_1, growth_rate=g_1, population_id=0),
-        msprime.PopulationParametersChange(time=t_2, growth_rate=g_2, population_id=0),
-        msprime.PopulationParametersChange(time=t_3, growth_rate=g_3, population_id=0),
-        msprime.PopulationParametersChange(time=t_4, growth_rate=g_4, population_id=0),
-        msprime.PopulationParametersChange(time=t_5, growth_rate=g_5, population_id=0),
-        msprime.PopulationParametersChange(
-            time=t_ancient, initial_size=n_ancient, growth_rate=0, population_id=0
-        ),
-    ]
-
-    return stdpopsim.DemographicModel(
-        id=id,
-        description=description,
-        long_description=long_description,
-        populations=populations,
-        citations=citations,
-        generation_time=generation_time,
-        population_configurations=population_configurations,
-        demographic_events=demographic_events,
-    )
-
-
-_species.add_demographic_model(_zigzag())
-
-
-def _kamm_ancient_eurasia():
-    id = "AncientEurasia_9K19"
-    description = "Multi-population model of ancient Eurasia"
-    long_description = """
-        This is the best-fitting model of a history of
-        multiple ancient and present-day human populations
-        sampled across Eurasia over the past 120,000 years.
-        The fitting was performed using momi2 (Kamm et al. 2019),
-        which uses the multi-population site-frequency spectrum
-        as input data. The model includes a ghost admixture event
-        (from unsampled basal Eurasians into early European
-        farmers), and two admixture events where the source is
-        approximately well-known (from Neanderthals into
-        Non-Africans and from Western European hunter-gatherers
-        into modern Sardinians. There are three present-day
-        populations: Sardinians, Han Chinese and African Mbuti.
-        Additionally, there are several ancient samples
-        obtained from fossils dated at different times in
-        the past: the Altai Neanderthal (Prufer et al. 2014),
-        a Mesolithic hunter-gatherer (Lazaridis et al. 2014),
-        a Neolithic early European sample (Lazaridis et al. 2014),
-        and two Palaeolithic modern humans from Siberia - MA1
-        (Raghavan et al. 2014) and  Ust'Ishim (Fu et al. 2014).
-        All the ancient samples are represented by a single diploid
-        genome.
-    """
-    # Sampling times are assuming 25 years per generation
-    populations = [
-        stdpopsim.Population(
-            id="Mbuti", description="Present-day African Mbuti", sampling_time=0
-        ),
-        # LBK: 8,000 years ago
-        stdpopsim.Population(
-            id="LBK", description="Early European farmer (EEF)", sampling_time=320
-        ),
-        stdpopsim.Population(
-            id="Sardinian", description="Present-day Sardinian", sampling_time=0
-        ),
-        # Loschbour: 7,500 years ago
-        stdpopsim.Population(
-            id="Loschbour",
-            description="Western hunter-gatherer (WHG)",
-            sampling_time=300,
-        ),
-        # MA1: 24,000 years ago
-        stdpopsim.Population(
-            id="MA1", description="Upper Palaeolithic MAl'ta culture", sampling_time=960
-        ),
-        stdpopsim.Population(
-            id="Han", description="Present-day Han Chinese", sampling_time=0
-        ),
-        # Ust Ishim: 45,000 years ago
-        stdpopsim.Population(
-            id="UstIshim",
-            description="early Palaeolithic Ust'-Ishim",
-            sampling_time=1800,
-        ),
-        # Altai Neanderthal: 50,000 years ago
-        stdpopsim.Population(
-            id="Neanderthal",
-            description="Altai Neanderthal from Siberia",
-            sampling_time=2000,
-        ),
-        stdpopsim.Population(
-            id="BasalEurasian", description="Basal Eurasians", sampling_time=None
-        ),
-    ]
-    citations = [
-        stdpopsim.Citation(
-            author="Kamm et al.",
-            year=2019,
-            doi="https://doi.org/10.1080/01621459.2019.1635482",
-            reasons={stdpopsim.CiteReason.DEM_MODEL},
-        )
-    ]
-
-    # Times are provided in years, so we convert into generations.
-    generation_time = 25
-    # Mutation_rate in Kamm et al. = 1.22e-8
-    # Effective population sizes
-    N_Losch = 1920
-    N_Mbu = 17300
-    N_Mbu_Losch = 29100
-    N_Han = 6300
-    N_Han_Losch = 2340
-    N_Nean_Losch = 18200
-    N_Nean = 86.9
-    N_LBK = 75.7
-    N_Sard = 15000
-    N_Sard_LBK = 12000
-    # Table A.1 has Altai at 50,000 years ago
-    t_NeaPopSizeChange = 50000 / generation_time
-    # Unknown but suspected parameters...
-    N_Basal = N_Losch
-    N_MA1 = N_Losch
-    N_Ust = N_Losch
-    # Population split times
-    t_Mbu_Losch = 95800 / generation_time
-    t_Han_Losch = 50400 / generation_time
-    t_Ust_Losch = 51500 / generation_time
-    t_Nean_Losch = 696000 / generation_time
-    t_MA1_Losch = 44900 / generation_time
-    t_LBK_Losch = 37700 / generation_time
-    t_Basal_Losch = 79800 / generation_time
-    t_Sard_LBK = 7690 / generation_time
-    # Given that we're using best model estimate,
-    # ghost WHG is directly descended from Loschbour,
-    # so this parameter is not used
-    # t_GhostWHG_Losch = 1560 / generation_time
-    # Admixture times
-    t_Nean_to_Eurasian = 56800 / generation_time
-    t_Basal_to_EEF = 33700 / generation_time
-    t_GhostWHG_to_Sard = 1230 / generation_time
-    t_NeanGrowth = t_Mbu_Losch - t_NeaPopSizeChange
-    logdiffNeanGrowth = math.log(N_Nean / N_Nean_Losch)
-    r_NeanGrowth = logdiffNeanGrowth / t_NeanGrowth
-    p_Nean_to_Eurasian = 0.0296
-    p_Basal_to_EEF = 0.0936
-    p_GhostWHG_to_Sard = 0.0317
-    # Population IDs: Mbuti = 0; LBK = 1;
-    # Sardinian = 2; Loschbour = 3; MA1 = 4;
-    # Han = 5; Ust Ishim = 6; Neanderthal = 7;
-    # Basal Eurasian = 8
-    population_configurations = [
-        msprime.PopulationConfiguration(
-            initial_size=N_Mbu, metadata=populations[0].asdict()
-        ),
-        msprime.PopulationConfiguration(
-            initial_size=N_LBK, metadata=populations[1].asdict()
-        ),
-        msprime.PopulationConfiguration(
-            initial_size=N_Sard, metadata=populations[2].asdict()
-        ),
-        msprime.PopulationConfiguration(
-            initial_size=N_Losch, metadata=populations[3].asdict()
-        ),
-        msprime.PopulationConfiguration(
-            initial_size=N_MA1, metadata=populations[4].asdict()
-        ),
-        msprime.PopulationConfiguration(
-            initial_size=N_Han, metadata=populations[5].asdict()
-        ),
-        msprime.PopulationConfiguration(
-            initial_size=N_Ust, metadata=populations[6].asdict()
-        ),
-        msprime.PopulationConfiguration(
-            initial_size=N_Nean, metadata=populations[7].asdict()
-        ),
-        msprime.PopulationConfiguration(
-            initial_size=N_Basal, metadata=populations[8].asdict()
-        ),
-    ]
-    demographic_events = [
-        # Sardinian receives admixture from Loschbour / WHG
-        msprime.MassMigration(
-            time=t_GhostWHG_to_Sard,
-            source=2,
-            destination=3,
-            proportion=p_GhostWHG_to_Sard,
-        ),
-        # Sardinian merges into LBK / EEF
-        # Now pop 1: Sardinian-LBK ancestral pop
-        msprime.MassMigration(time=t_Sard_LBK, source=2, destination=1, proportion=1.0),
-        # Sardinian-LBK ancestral pop size change
-        msprime.PopulationParametersChange(
-            time=t_Sard_LBK, initial_size=N_Sard_LBK, population_id=1
-        ),
-        # LBK / EEF receives admixture from Basal Eurasians
-        msprime.MassMigration(
-            time=t_Basal_to_EEF, source=1, destination=8, proportion=p_Basal_to_EEF
-        ),
-        # LBK / EEF merges into Loschbour
-        msprime.MassMigration(
-            time=t_LBK_Losch, source=1, destination=3, proportion=1.0
-        ),
-        # MA1 merges into Loschbour
-        msprime.MassMigration(
-            time=t_MA1_Losch, source=4, destination=3, proportion=1.0
-        ),
-        # Neanderthal start change in population size
-        msprime.PopulationParametersChange(
-            time=t_NeaPopSizeChange,
-            initial_size=N_Nean,
-            growth_rate=r_NeanGrowth,
-            population_id=7,
-        ),
-        # Han merges into Loschbour
-        msprime.MassMigration(
-            time=t_Han_Losch, source=5, destination=3, proportion=1.0
-        ),
-        # Change in population size in Han-Losch ancestral pop
-        msprime.PopulationParametersChange(
-            time=t_Han_Losch, initial_size=N_Han_Losch, population_id=3
-        ),
-        # UstIshim merges into Loschbour
-        msprime.MassMigration(
-            time=t_Ust_Losch, source=6, destination=3, proportion=1.0
-        ),
-        # Loschbour / Non-Africans receive admixture from Neanderthals
-        msprime.MassMigration(
-            time=t_Nean_to_Eurasian,
-            source=3,
-            destination=7,
-            proportion=p_Nean_to_Eurasian,
-        ),
-        # Basal Eurasians merge into Loschbour / Non-Africans
-        msprime.MassMigration(
-            time=t_Basal_Losch, source=8, destination=3, proportion=1.0
-        ),
-        # Mbuti merge into Loschbour / Non-Africans
-        msprime.MassMigration(
-            time=t_Mbu_Losch, source=0, destination=3, proportion=1.0
-        ),
-        # Change in population size in Mbuti-Losch ancestral pop
-        msprime.PopulationParametersChange(
-            time=t_Mbu_Losch, initial_size=N_Mbu_Losch, population_id=3
-        ),
-        # Change in population size in Neanderthal, growth rate 0
-        msprime.PopulationParametersChange(
-            time=t_Mbu_Losch, initial_size=N_Nean_Losch, growth_rate=0, population_id=7
-        ),
-        # Neanderthal merge into Loschbour / modern humans
-        msprime.MassMigration(
-            time=t_Nean_Losch, source=7, destination=3, proportion=1.0
-        ),
-        # Ancestral hominin population size change
-        msprime.PopulationParametersChange(
-            time=t_Nean_Losch, initial_size=N_Nean_Losch, population_id=3
-        ),
-    ]
-
-    return stdpopsim.DemographicModel(
-        id=id,
-        description=description,
-        long_description=long_description,
-        populations=populations,
-        citations=citations,
-        generation_time=generation_time,
-        population_configurations=population_configurations,
-        demographic_events=demographic_events,
-    )
-
-
-_species.add_demographic_model(_kamm_ancient_eurasia())
-
-
-def _papuans_10j19():
-    id = "PapuansOutOfAfrica_10J19"
-    description = "Out-of-Africa with archaic admixture into Papuans"
-    long_description = """
-        A ten population model of out-of-Africa, including two pulses of
-        Denisovan admixture into Papuans, and several pulses of Neandertal
-        admixture into non-Africans.
-        Most parameters are from Jacobs et al. (2019), Table S5 and Figure S5.
-        This model is an extension of one from Malaspinas et al. (2016), thus
-        some parameters are inherited from there.
-        """
-
-    # sampling times
-    T_DenA = 2058
-    T_NeaA = 2612
-
-    populations = [
-        # humans
-        _yri_population,
-        _ceu_population,
-        _chb_population,
-        stdpopsim.Population("Papuan", "Papuans from Indonesia and New Guinea"),
-        # archaics
-        stdpopsim.Population(
-            "DenA", "Altai Denisovan (sampling) lineage", sampling_time=T_DenA
-        ),
-        stdpopsim.Population(
-            "NeaA", "Altai Neandertal (sampling) lineage", sampling_time=T_NeaA
-        ),
-        stdpopsim.Population(
-            "Den1", "Denisovan D1 (introgressing) lineage", sampling_time=None
-        ),
-        stdpopsim.Population(
-            "Den2", "Denisovan D2 (introgressing) lineage", sampling_time=None
-        ),
-        stdpopsim.Population(
-            "Nea1", "Neandertal N1 (introgressing) lineage", sampling_time=None
-        ),
-        stdpopsim.Population("Ghost", "Out-of-Africa lineage", sampling_time=None),
-    ]
-    pop = {p.id: i for i, p in enumerate(populations)}
-
-    citations = [
-        stdpopsim.Citation(
-            author="Jacobs et al.",
-            year=2019,
-            doi="https://doi.org/10.1016/j.cell.2019.02.035",
-            reasons={stdpopsim.CiteReason.DEM_MODEL},
-        ),
-        stdpopsim.Citation(
-            author="Malaspinas et al.",
-            year=2016,
-            doi="https://doi.org/10.1038/nature18299",
-            reasons={stdpopsim.CiteReason.DEM_MODEL},
-        ),
-    ]
-
-    # Inherited from Malaspinas et al., which gives the following refs:
-    generation_time = 29  # Fenner 2005
-    # mutation_rate = 1.25e-8  # per gen per site, Scally & Durbin 2012
-
-    N_YRI = 48433
-    N_Ghost = 8516
-    N_CEU = 6962
-    N_CHB = 9025
-    N_Papuan = 8834
-
-    N_DenA = 5083
-    N_Den1 = 13249
-    N_Den2 = 13249
-    N_NeaA = 826
-    N_Nea1 = 13249
-
-    pop_meta = {p.id: p.asdict() for p in populations}
-    population_configurations = [
-        msprime.PopulationConfiguration(initial_size=N_YRI, metadata=pop_meta["YRI"]),
-        msprime.PopulationConfiguration(initial_size=N_CEU, metadata=pop_meta["CEU"]),
-        msprime.PopulationConfiguration(initial_size=N_CHB, metadata=pop_meta["CHB"]),
-        msprime.PopulationConfiguration(
-            initial_size=N_Papuan, metadata=pop_meta["Papuan"]
-        ),
-        msprime.PopulationConfiguration(initial_size=N_DenA, metadata=pop_meta["DenA"]),
-        msprime.PopulationConfiguration(initial_size=N_NeaA, metadata=pop_meta["NeaA"]),
-        msprime.PopulationConfiguration(initial_size=N_Den1, metadata=pop_meta["Den1"]),
-        msprime.PopulationConfiguration(initial_size=N_Den2, metadata=pop_meta["Den2"]),
-        msprime.PopulationConfiguration(initial_size=N_Nea1, metadata=pop_meta["Nea1"]),
-        msprime.PopulationConfiguration(
-            initial_size=N_Ghost, metadata=pop_meta["Ghost"]
-        ),
-    ]
-    assert len(populations) == len(population_configurations)
-
-    # initial migrations
-    m_Ghost_Afr = 1.79e-4
-    m_Ghost_EU = 4.42e-4
-    m_EU_AS = 3.14e-5
-    m_AS_Papuan = 5.72e-5
-    # older migrations
-    m_Eurasian_Papuan = 5.72e-4
-    m_Eurasian_Ghost = 4.42e-4
-
-    migration_matrix = [[0] * len(populations) for _ in range(len(populations))]
-    migration_matrix[pop["Ghost"]][pop["YRI"]] = m_Ghost_Afr
-    migration_matrix[pop["YRI"]][pop["Ghost"]] = m_Ghost_Afr
-    migration_matrix[pop["Ghost"]][pop["CEU"]] = m_Ghost_EU
-    migration_matrix[pop["CEU"]][pop["Ghost"]] = m_Ghost_EU
-    migration_matrix[pop["CEU"]][pop["CHB"]] = m_EU_AS
-    migration_matrix[pop["CHB"]][pop["CEU"]] = m_EU_AS
-    migration_matrix[pop["CHB"]][pop["Papuan"]] = m_AS_Papuan
-    migration_matrix[pop["Papuan"]][pop["CHB"]] = m_AS_Papuan
-
-    # splits
-    T_EU_AS_split = 1293
-    T_Eurasian_Ghost_split = 1758
-    T_Papuan_Ghost_split = 1784
-    T_Ghost_Afr_split = 2218
-    T_NeaA_Nea1_split = 3375
-    T_DenA_Den1_split = 9750
-    T_DenA_Den2_split = 12500
-    T_Den_Nea_split = 15090
-    T_Afr_Archaic_split = 20225
-
-    # bottlenecks
-    Tb_Eurasia = 1659
-    Tb_Papua = 1685
-    Tb_Ghost = 2119
-    Nb_Eurasia = 2231
-    Nb_Papua = 243
-    Nb_Ghost = 1394
-
-    # internal branches
-    N_EU_AS = 12971
-    N_Ghost_Afr = 41563
-    N_NeaA_Nea1 = 13249
-    N_Den_Anc = 100  # S10.i p. 31/45
-    N_DenA_Den1 = N_Den_Anc
-    N_DenA_Den2 = N_Den_Anc
-    N_Den_Nea = 13249
-    N_Afr_Archaic = 32671
-
-    # admixture pulses
-    m_Den_Papuan = 0.04
-    p1 = 0.55  # S10.i p. 31/45
-    m_Den1_Papuan = p1 * m_Den_Papuan
-    m_Den2_Papuan = (1 - p1) * m_Den_Papuan
-    m_Nea1_Ghost = 0.024
-    m_Nea1_Eurasian = 0.011
-    m_Nea1_Papuan = 0.002
-    m_Nea1_AS = 0.002
-
-    T_Nea1_Ghost_mig = 1853
-    T_Nea1_Eurasian_mig = 1566
-    T_Nea1_Papuan_mig = 1412
-    T_Nea1_AS_mig = 883
-    # Fig. 4B, and S10.h p. 30/45
-    T_Den1_Papuan_mig = 29.8e3 / generation_time
-    T_Den2_Papuan_mig = 45.7e3 / generation_time
-
-    demographic_events = [
-        # human lineage splits
-        msprime.MassMigration(
-            time=T_EU_AS_split, source=pop["CEU"], destination=pop["CHB"]
-        ),
-        msprime.PopulationParametersChange(
-            time=T_EU_AS_split, initial_size=N_EU_AS, population_id=pop["CHB"]
-        ),
-        msprime.MassMigration(
-            time=T_Eurasian_Ghost_split, source=pop["CHB"], destination=pop["Ghost"]
-        ),
-        msprime.MassMigration(
-            time=T_Papuan_Ghost_split, source=pop["Papuan"], destination=pop["Ghost"]
-        ),
-        msprime.MassMigration(
-            time=T_Ghost_Afr_split, source=pop["Ghost"], destination=pop["YRI"]
-        ),
-        msprime.PopulationParametersChange(
-            time=T_Ghost_Afr_split, initial_size=N_Ghost_Afr, population_id=pop["YRI"]
-        ),
-        # archaic lineage splits
-        msprime.MassMigration(
-            time=T_DenA_Den1_split, source=pop["Den1"], destination=pop["DenA"]
-        ),
-        msprime.PopulationParametersChange(
-            time=T_DenA_Den1_split, initial_size=N_DenA_Den1, population_id=pop["DenA"]
-        ),
-        msprime.MassMigration(
-            time=T_DenA_Den2_split, source=pop["Den2"], destination=pop["DenA"]
-        ),
-        msprime.PopulationParametersChange(
-            time=T_DenA_Den2_split, initial_size=N_DenA_Den2, population_id=pop["DenA"]
-        ),
-        msprime.MassMigration(
-            time=T_NeaA_Nea1_split, source=pop["Nea1"], destination=pop["NeaA"]
-        ),
-        msprime.PopulationParametersChange(
-            time=T_NeaA_Nea1_split, initial_size=N_NeaA_Nea1, population_id=pop["NeaA"]
-        ),
-        msprime.MassMigration(
-            time=T_Den_Nea_split, source=pop["NeaA"], destination=pop["DenA"]
-        ),
-        msprime.PopulationParametersChange(
-            time=T_Den_Nea_split, initial_size=N_Den_Nea, population_id=pop["DenA"]
-        ),
-        msprime.MassMigration(
-            time=T_Afr_Archaic_split, source=pop["DenA"], destination=pop["YRI"]
-        ),
-        msprime.PopulationParametersChange(
-            time=T_Afr_Archaic_split,
-            initial_size=N_Afr_Archaic,
-            population_id=pop["YRI"],
-        ),
-        # bottlenecks
-        msprime.PopulationParametersChange(
-            time=Tb_Eurasia, initial_size=Nb_Eurasia, population_id=pop["CHB"]
-        ),
-        msprime.PopulationParametersChange(
-            time=Tb_Papua, initial_size=Nb_Papua, population_id=pop["Papuan"]
-        ),
-        msprime.PopulationParametersChange(
-            time=Tb_Ghost, initial_size=Nb_Ghost, population_id=pop["Ghost"]
-        ),
-        # migration changes
-        msprime.MigrationRateChange(
-            time=T_EU_AS_split, rate=0, matrix_index=(pop["CHB"], pop["CEU"])
-        ),
-        msprime.MigrationRateChange(
-            time=T_EU_AS_split, rate=0, matrix_index=(pop["CEU"], pop["CHB"])
-        ),
-        msprime.MigrationRateChange(
-            time=T_EU_AS_split, rate=0, matrix_index=(pop["Papuan"], pop["CHB"])
-        ),
-        msprime.MigrationRateChange(
-            time=T_EU_AS_split, rate=0, matrix_index=(pop["CHB"], pop["Papuan"])
-        ),
-        msprime.MigrationRateChange(
-            time=T_EU_AS_split, rate=0, matrix_index=(pop["Ghost"], pop["CEU"])
-        ),
-        msprime.MigrationRateChange(
-            time=T_EU_AS_split, rate=0, matrix_index=(pop["CEU"], pop["Ghost"])
-        ),
-        msprime.MigrationRateChange(
-            time=T_EU_AS_split,
-            rate=m_Eurasian_Papuan,
-            matrix_index=(pop["CHB"], pop["Papuan"]),
-        ),
-        msprime.MigrationRateChange(
-            time=T_EU_AS_split,
-            rate=m_Eurasian_Papuan,
-            matrix_index=(pop["Papuan"], pop["CHB"]),
-        ),
-        msprime.MigrationRateChange(
-            time=T_EU_AS_split,
-            rate=m_Eurasian_Ghost,
-            matrix_index=(pop["CHB"], pop["Ghost"]),
-        ),
-        msprime.MigrationRateChange(
-            time=T_EU_AS_split,
-            rate=m_Eurasian_Ghost,
-            matrix_index=(pop["Ghost"], pop["CHB"]),
-        ),
-        # all migrations off
-        msprime.MigrationRateChange(time=Tb_Eurasia, rate=0),
-        # admixture pulses
-        msprime.MassMigration(
-            time=T_Den1_Papuan_mig,
-            proportion=m_Den1_Papuan,
-            source=pop["Papuan"],
-            destination=pop["Den1"],
-        ),
-        msprime.MassMigration(
-            time=T_Den2_Papuan_mig,
-            proportion=m_Den2_Papuan,
-            source=pop["Papuan"],
-            destination=pop["Den2"],
-        ),
-        msprime.MassMigration(
-            time=T_Nea1_Ghost_mig,
-            proportion=m_Nea1_Ghost,
-            source=pop["Ghost"],
-            destination=pop["Nea1"],
-        ),
-        msprime.MassMigration(
-            time=T_Nea1_Eurasian_mig,
-            proportion=m_Nea1_Eurasian,
-            source=pop["CHB"],
-            destination=pop["Nea1"],
-        ),
-        msprime.MassMigration(
-            time=T_Nea1_Papuan_mig,
-            proportion=m_Nea1_Papuan,
-            source=pop["Papuan"],
-            destination=pop["Nea1"],
-        ),
-        msprime.MassMigration(
-            time=T_Nea1_AS_mig,
-            proportion=m_Nea1_AS,
-            source=pop["CHB"],
-            destination=pop["Nea1"],
-        ),
-    ]
-
-    demographic_events.sort(key=lambda x: x.time)
-
-    return stdpopsim.DemographicModel(
-        id=id,
-        description=description,
-        long_description=long_description,
-        populations=populations,
-        citations=citations,
-        generation_time=generation_time,
-        population_configurations=population_configurations,
-        migration_matrix=migration_matrix,
-        demographic_events=demographic_events,
-    )
-
-
-_species.add_demographic_model(_papuans_10j19())
-
-
-def _AJ():
-    id = "AshkSub_7G19"
-    description = "Ashkenazi Jewish with substructure and European admixture"
-    long_description = """
-    This was the best fit model of Ashkenazi Jewish demographic history from
-    Gladstein and Hammer 2019, shown in Figure 1, labeled "Substructure Model".
-    Model choice and parameter estimation were performed with Approximate
-    Bayesian Computation. Parameter values are based on the mode from ABC found
-    in Table S3 of Gladstein and Hammer 2019. In this model, the ancestors of
-    Europeans and Middle Eastern populations diverge. Non-Ashkenazi Jewish
-    populations then diverge from the Middle Eastern population. The Ashkenazi
-    Jews then diverge from the other Jewish populations and experience a
-    substantial reduction in population size and a single pulse of gene flow
-    from Europeans (corresponding to their arrival in Europe). After the gene
-    flow from Europeans to the Ashkenazi Jews, the Ashkenazi Jews split into
-    two groups, the Western and Eastern. Finally, the Western Ashkenazi Jews
-    experience moderate instantaneous population size increase, and the
-    Eastern experience a massive population size increase. In addition to the
-    demographic model Gladstein and Hammer 2019 also incorporated an SNP array
-    ascertainment scheme into the simulation. This demographic model does not
-    include the SNP array ascertainment scheme. It should be noted that
-    Gladstein and Hammer 2019 simulated with a mutation rate of 2.5e-8.
-    """
-    populations = [
-        _yri_population,
-        _chb_population,
-        _ceu_population,
-        stdpopsim.Population(id="ME", description="Middle Eastern"),
-        stdpopsim.Population(id="J", description="non-Ashkenazi Jewish"),
-        stdpopsim.Population(id="WAJ", description="Western Ashkenazi Jewish"),
-        stdpopsim.Population(id="EAJ", description="Eastern Ashkenazi Jewish"),
-    ]
-    citations = [
-        stdpopsim.Citation(
-            author="Gladstein and Hammer",
-            year=2019,
-            doi="https://doi.org/10.1093/molbev/msz047",
-            reasons={stdpopsim.CiteReason.DEM_MODEL},
-        )
-    ]
-
-    generation_time = 25
-
-    # parameter value definitions based on mode from ABC
-    # found in Table S3 of Gladstein and Hammer 2019
-
-    # effective population sizes
-    NANC = 7300  # not inferred. Value taken from Gutenkunst et al. 2009
-    NYRI = 10 ** 4.26
-    NCHB = 10 ** 3.61
-    NCEU = 10 ** 4.52
-    NM = 10 ** 5.64
-    NJ = 10 ** 5.55
-    NAg = 10 ** 3.04
-    NWA = 10 ** 3.82
-    NEA = 10 ** 6.29
-
-    # admixture proportion from European to Ashkenazi Jews.
-    m = 0.17
-
-    # Times in generations
-    Tgrowth = 8800  # not inferred. Value taken from Gutenkunst et al. 2009
-    TAF = 2105
-    Teu_as = 850
-    TEM = 481
-    TMJ = 211
-    TA = 29
-    Tm = 28  # not inferred. For simplicity put at one generation before TA
-    TAEW = 14
-    TAg = 13  # not inferred. For simplicity put at one generation before TAEW
-
-    YRI, CHB, CEU, M, J, WA, EA = 0, 1, 2, 3, 4, 5, 6
-
-    return stdpopsim.DemographicModel(
-        id=id,
-        description=description,
-        long_description=long_description,
-        populations=populations,
-        citations=citations,
-        generation_time=generation_time,
-        population_configurations=[
-            msprime.PopulationConfiguration(
-                initial_size=NYRI, metadata=populations[0].asdict()
-            ),
-            msprime.PopulationConfiguration(
-                initial_size=NCHB, metadata=populations[1].asdict()
-            ),
-            msprime.PopulationConfiguration(
-                initial_size=NCEU, metadata=populations[2].asdict()
-            ),
-            msprime.PopulationConfiguration(
-                initial_size=NM, metadata=populations[3].asdict()
-            ),
-            msprime.PopulationConfiguration(
-                initial_size=NJ, metadata=populations[4].asdict()
-            ),
-            msprime.PopulationConfiguration(
-                initial_size=NWA, metadata=populations[5].asdict()
-            ),
-            msprime.PopulationConfiguration(
-                initial_size=NEA, metadata=populations[6].asdict()
-            ),
-        ],
-        demographic_events=[
-            # instantaneous growth in EA and WA at the same time
-            msprime.PopulationParametersChange(
-                time=TAg, initial_size=NAg, population_id=WA
-            ),
-            msprime.PopulationParametersChange(
-                time=TAg, initial_size=NAg, population_id=EA
-            ),
-            # EA splits from WA
-            msprime.MassMigration(time=TAEW, source=EA, destination=WA, proportion=1.0),
-            # E geneflow into WA ancestor (forward in time)
-            msprime.MassMigration(time=Tm, source=WA, destination=CEU, proportion=m),
-            # WA ancestor splits from J
-            msprime.MassMigration(time=TA, source=WA, destination=J, proportion=1.0),
-            # J splits from M split
-            msprime.MassMigration(time=TMJ, source=J, destination=M, proportion=1.0),
-            # M splits from CEU
-            msprime.MassMigration(time=TEM, source=M, destination=CEU, proportion=1.0),
-            # CEU splits from CHB
-            msprime.MassMigration(
-                time=Teu_as, source=CEU, destination=CHB, proportion=1.0
-            ),
-            # CHB splits from YRI
-            msprime.MassMigration(
-                time=TAF, source=CHB, destination=YRI, proportion=1.0
-            ),
-            # Instantaneous growth in YRI
-            msprime.PopulationParametersChange(
-                time=Tgrowth, initial_size=NANC, population_id=YRI
-            ),
-        ],
-    )
-
-
-_species.add_demographic_model(_AJ())
-
-
-def _ooa_4pop():
-    id = "OutOfAfrica_4J17"
-    description = "4 population out of Africa"
-    long_description = """
-        Demographic model for a four population out-of-Africa history,
-        taken from Jouganous et al. (2017). Parameter values were taken
-        from table 4 in the main text. This model was fit based on joint
-        allele frequecy spectrum (AFS) data from 1000 Genomes exomes
-        from the YRI, CEU, CHB, and JPT poulation samples. The demography
-        follows the previous three-populations out-of-Africa models
-        with an additional population split in Asia leading to the
-        Japanese (JPT) population. Parameter values were estimated with
-        the program Moments assuming a mutation rate of 1.44e-8 and a
-        generation time of 29 years.
-    """
-    populations = [
-        _yri_population,
-        _ceu_population,
-        _chb_population,
-        stdpopsim.Population(
-            id="JPT",
-            description="1000 Genomes JPT (Japanese in Tokyo, Japan)",
-            sampling_time=0,
-        ),
-    ]
-    citations = [
-        stdpopsim.Citation(
-            author="Jouganous et al.",
-            year="2017",
-            doi="https://doi.org/10.1534/genetics.117.200493",
-            reasons={stdpopsim.CiteReason.DEM_MODEL},
-        )
-    ]
-
-    generation_time = 29
-
-    # Parameter values from Table 4 (in bold)
-
-    # Population sizes
-    N_A = 11293  # ancestral YRI population size
-    N_AF = 23721  # YRI population size
-    N_B = 2831  # B population size (out-of-Africa)
-    N_EU0 = 2512  # CEU population size (intial)
-    N_AS0 = 1019  # Asian population size (inital)
-    N_JP0 = 4384  # Japanese population size (inital)
-
-    # growth rates
-    r_EU = 0.0016  # European growth rate
-    r_JP = 0.0129  # Japanese growth rate
-    r_AS = 0.0026  # Asian growth rate
-
-    # migration rates
-    m_AF_B = 16.8e-5  # YRI <--> B (out-of-Africa)
-    m_AF_EU = 1.14e-5  # YRI <--> CEU
-    m_AF_AS = 0.56e-5  # YRI <--> CHB
-    m_EU_AS = 4.75e-5  # European <--> CHB
-    m_CH_JP = 3.3e-5  # CHBHan <--> JPT
-
-    # epochs (specified in years, converted to generations)
-    T_AF = int(357000 / generation_time)  # time ancestral to African
-    T_B = int(119000 / generation_time)  # time of Out of Africa split
-    T_EU_AS = int(46000 / generation_time)  # time of Europe/Asia split
-    T_CH_JP = int(9000 / generation_time)  # time of Han/Japanese split
-
-    # also
-    r_AF = 0  # YRI and B growth rate
-
-    # A and YRI
-    pop0 = msprime.PopulationConfiguration(
-        initial_size=N_AF * math.exp(r_AF * T_AF),
-        growth_rate=r_AF,
-        metadata=populations[0].asdict(),
-    )
-
-    # B and CEU
-    pop1 = msprime.PopulationConfiguration(
-        initial_size=N_EU0 * math.exp(r_EU * T_EU_AS),
-        growth_rate=r_EU,
-        metadata=populations[1].asdict(),
-    )
-
-    # Asian and CHB
-    pop2 = msprime.PopulationConfiguration(
-        initial_size=N_AS0 * math.exp(r_AS * T_EU_AS),
-        growth_rate=r_AS,
-        metadata=populations[2].asdict(),
-    )
-
-    # JPT
-    pop3 = msprime.PopulationConfiguration(
-        initial_size=N_JP0 * math.exp(r_JP * T_CH_JP),
-        growth_rate=r_JP,
-        metadata=populations[3].asdict(),
-    )
-
-    population_configurations = [pop0, pop1, pop2, pop3]
-
-    migration_matrix = [
-        [0, m_AF_EU, m_AF_AS, 0],
-        [m_AF_EU, 0, m_EU_AS, 0],
-        [m_AF_AS, m_EU_AS, 0, m_CH_JP],
-        [0, 0, m_CH_JP, 0],
-    ]
-
-    # Demographic events from most recent -> oldest
-
-    # JPT splits from CHB
-    JPT_event = [
-        # seg migration rate JPT <--> CHB to zero
-        msprime.MigrationRateChange(time=T_CH_JP, rate=0.0, matrix_index=(2, 3)),
-        msprime.MigrationRateChange(time=T_CH_JP, rate=0.0, matrix_index=(3, 2)),
-        # JPT splits from CHB
-        msprime.MassMigration(time=T_CH_JP, source=3, destination=2, proportion=1.0),
-    ]
-
-    # CEU and CHB merge into B
-    AS_event = [
-        # set migration rates to zero
-        msprime.MigrationRateChange(time=T_EU_AS, rate=0.0),
-        # CHB splits from CEU
-        msprime.MassMigration(time=T_EU_AS, source=2, destination=1, proportion=1.0),
-        # change size of B to N_B
-        # set growth rates in B to zero
-        msprime.PopulationParametersChange(
-            time=T_EU_AS, initial_size=N_B, growth_rate=0.0, population_id=1
-        ),
-        # set up migration between YRI and B
-        msprime.MigrationRateChange(time=T_EU_AS, rate=m_AF_B, matrix_index=(0, 1)),
-        msprime.MigrationRateChange(time=T_EU_AS, rate=m_AF_B, matrix_index=(1, 0)),
-    ]
-
-    # Out-of-Africa event
-    EU_event = [
-        # B splits from YRI
-        msprime.MigrationRateChange(time=T_B, rate=0.0),
-        msprime.MassMigration(time=T_B, source=1, destination=0, proportion=1.0),
-    ]
-
-    # Change size of YRI/African population
-    AF_event = [
-        msprime.PopulationParametersChange(time=T_AF, initial_size=N_A, population_id=0)
-    ]
-
-    demographic_events = JPT_event + AS_event + EU_event + AF_event
-
-    return stdpopsim.DemographicModel(
-        id=id,
-        description=description,
-        long_description=long_description,
-        populations=populations,
-        citations=citations,
-        generation_time=generation_time,
-        population_configurations=population_configurations,
-        migration_matrix=migration_matrix,
-        demographic_events=demographic_events,
-    )
-
-
-_species.add_demographic_model(_ooa_4pop())
+from . import species  # noqa 501
+from . import genetic_maps  # noqa 501
+from . import demographic_models  # noqa 501
+from . import annotations  # noqa 501
diff --git a/stdpopsim/catalog/HomSap/annotations.py b/stdpopsim/catalog/HomSap/annotations.py
new file mode 100644
index 000000000..27f98117a
--- /dev/null
+++ b/stdpopsim/catalog/HomSap/annotations.py
@@ -0,0 +1,26 @@
+import stdpopsim
+
+_species = stdpopsim.get_species("HomSap")
+
+_an = stdpopsim.Annotation(
+    species=_species,
+    id="Ensembl_GRCh38_gff3",
+    description="Ensembl GFF3 annotations on GRCh38",
+    url=(
+        "ftp://ftp.ensembl.org/pub/release-101/"
+        "gff3/homo_sapiens/Homo_sapiens.GRCh38.101.gff3.gz"
+    ),
+    zarr_url=(
+        "https://stdpopsim.s3-us-west-2.amazonaws.com/annotations/HomSap.GRCh38.zip"
+    ),
+    zarr_sha256="929c52dc5e5172d4a96a776d066f014f8f207f1477a21f7b88626f81a7142d0b",
+    citations=[
+        stdpopsim.Citation(
+            year=2018,
+            author="Hunt et al",
+            doi="https://doi.org/10.1093/database/bay119",
+            reasons={stdpopsim.CiteReason.ANNOTATION},
+        )
+    ],
+)
+_species.add_annotations(_an)
diff --git a/stdpopsim/catalog/HomSap/demographic_models.py b/stdpopsim/catalog/HomSap/demographic_models.py
new file mode 100644
index 000000000..4e5906853
--- /dev/null
+++ b/stdpopsim/catalog/HomSap/demographic_models.py
@@ -0,0 +1,1516 @@
+import math
+import msprime
+
+import stdpopsim
+
+_species = stdpopsim.get_species("HomSap")
+
+###########################################################
+#
+# Demographic models
+#
+###########################################################
+
+# population definitions that are reused.
+_yri_population = stdpopsim.Population(
+    id="YRI", description="1000 Genomes YRI (Yorubans)"
+)
+_ceu_population = stdpopsim.Population(
+    id="CEU",
+    description=(
+        "1000 Genomes CEU (Utah Residents (CEPH) with Northern and "
+        "Western European Ancestry"
+    ),
+)
+_chb_population = stdpopsim.Population(
+    id="CHB", description="1000 Genomes CHB (Han Chinese in Beijing, China)"
+)
+
+
+_tennessen_et_al = stdpopsim.Citation(
+    author="Tennessen et al.",
+    year=2012,
+    doi="https://doi.org/10.1126/science.1219240",
+    reasons={stdpopsim.CiteReason.DEM_MODEL},
+)
+
+
+def _ooa_3():
+    id = "OutOfAfrica_3G09"
+    description = "Three population out-of-Africa"
+    long_description = """
+        The three population Out-of-Africa model from Gutenkunst et al. 2009.
+        It describes the ancestral human population in Africa, the out of Africa
+        event, and the subsequent European-Asian population split.
+        Model parameters are the maximum likelihood values of the
+        various parameters given in Table 1 of Gutenkunst et al.
+    """
+    populations = [_yri_population, _ceu_population, _chb_population]
+
+    citations = [
+        stdpopsim.Citation(
+            author="Gutenkunst et al.",
+            year=2009,
+            doi="https://doi.org/10.1371/journal.pgen.1000695",
+            reasons={stdpopsim.CiteReason.DEM_MODEL},
+        )
+    ]
+
+    generation_time = 25
+
+    # First we set out the maximum likelihood values of the various parameters
+    # given in Table 1.
+    N_A = 7300
+    N_B = 2100
+    N_AF = 12300
+    N_EU0 = 1000
+    N_AS0 = 510
+    # Times are provided in years, so we convert into generations.
+
+    T_AF = 220e3 / generation_time
+    T_B = 140e3 / generation_time
+    T_EU_AS = 21.2e3 / generation_time
+    # We need to work out the starting (diploid) population sizes based on
+    # the growth rates provided for these two populations
+    r_EU = 0.004
+    r_AS = 0.0055
+    N_EU = N_EU0 / math.exp(-r_EU * T_EU_AS)
+    N_AS = N_AS0 / math.exp(-r_AS * T_EU_AS)
+    # Migration rates during the various epochs.
+    m_AF_B = 25e-5
+    m_AF_EU = 3e-5
+    m_AF_AS = 1.9e-5
+    m_EU_AS = 9.6e-5
+
+    return stdpopsim.DemographicModel(
+        id=id,
+        description=description,
+        long_description=long_description,
+        populations=populations,
+        citations=citations,
+        generation_time=generation_time,
+        # Population IDs correspond to their indexes in the population
+        # configuration array. Therefore, we have 0=YRI, 1=CEU and 2=CHB
+        # initially.
+        population_configurations=[
+            msprime.PopulationConfiguration(
+                initial_size=N_AF, metadata=populations[0].asdict()
+            ),
+            msprime.PopulationConfiguration(
+                initial_size=N_EU, growth_rate=r_EU, metadata=populations[1].asdict()
+            ),
+            msprime.PopulationConfiguration(
+                initial_size=N_AS, growth_rate=r_AS, metadata=populations[2].asdict()
+            ),
+        ],
+        migration_matrix=[
+            [0, m_AF_EU, m_AF_AS],  # noqa
+            [m_AF_EU, 0, m_EU_AS],  # noqa
+            [m_AF_AS, m_EU_AS, 0],  # noqa
+        ],
+        demographic_events=[
+            # CEU and CHB merge into B with rate changes at T_EU_AS
+            msprime.MassMigration(
+                time=T_EU_AS, source=2, destination=1, proportion=1.0
+            ),
+            msprime.MigrationRateChange(time=T_EU_AS, rate=0),
+            msprime.MigrationRateChange(time=T_EU_AS, rate=m_AF_B, matrix_index=(0, 1)),
+            msprime.MigrationRateChange(time=T_EU_AS, rate=m_AF_B, matrix_index=(1, 0)),
+            msprime.PopulationParametersChange(
+                time=T_EU_AS, initial_size=N_B, growth_rate=0, population_id=1
+            ),
+            # Population B merges into YRI at T_B
+            msprime.MassMigration(time=T_B, source=1, destination=0, proportion=1.0),
+            msprime.MigrationRateChange(time=T_B, rate=0),
+            # Size changes to N_A at T_AF
+            msprime.PopulationParametersChange(
+                time=T_AF, initial_size=N_A, population_id=0
+            ),
+        ],
+    )
+
+
+_species.add_demographic_model(_ooa_3())
+
+
+def _ooa_2():
+    id = "OutOfAfrica_2T12"
+    description = "Two population out-of-Africa"
+    long_description = """
+        The model is derived from the Tennesen et al. analysis of the
+        jSFS from European Americans and African Americans.
+        It describes the ancestral human population in Africa, the out of Africa event,
+        and two distinct periods of subsequent European population growth over the past
+        23kya. Model parameters are taken from Fig. S5 in Fu et al.
+    """
+    populations = [
+        stdpopsim.Population(id="AFR", description="African Americans"),
+        stdpopsim.Population(id="EUR", description="European Americans"),
+    ]
+    citations = [
+        _tennessen_et_al,
+        stdpopsim.Citation(
+            author="Fu et al.",
+            year=2013,
+            doi="https://doi.org/10.1038/nature11690",
+            reasons={stdpopsim.CiteReason.DEM_MODEL},
+        ),
+    ]
+
+    generation_time = 25
+
+    T_AF = 148e3 / generation_time
+    T_OOA = 51e3 / generation_time
+    T_EU0 = 23e3 / generation_time
+    T_EG = 5115 / generation_time
+
+    # Growth rates
+    r_EU0 = 0.00307
+    r_EU = 0.0195
+    r_AF = 0.0166
+
+    # population sizes
+    N_A = 7310
+    N_AF1 = 14474
+    N_B = 1861
+    N_EU0 = 1032
+    N_EU1 = N_EU0 / math.exp(-r_EU0 * (T_EU0 - T_EG))
+
+    # migration rates
+    m_AF_B = 15e-5
+    m_AF_EU = 2.5e-5
+
+    # present Ne
+    N_EU = N_EU1 / math.exp(-r_EU * T_EG)
+    N_AF = N_AF1 / math.exp(-r_AF * T_EG)
+
+    return stdpopsim.DemographicModel(
+        id=id,
+        description=description,
+        long_description=long_description,
+        populations=populations,
+        citations=citations,
+        generation_time=generation_time,
+        population_configurations=[
+            msprime.PopulationConfiguration(
+                initial_size=N_AF, growth_rate=r_AF, metadata=populations[0].asdict()
+            ),
+            msprime.PopulationConfiguration(
+                initial_size=N_EU, growth_rate=r_EU, metadata=populations[1].asdict()
+            ),
+        ],
+        migration_matrix=[
+            [0, m_AF_EU],
+            [m_AF_EU, 0],
+        ],
+        demographic_events=[
+            msprime.MigrationRateChange(time=T_EG, rate=m_AF_EU, matrix_index=(0, 1)),
+            msprime.MigrationRateChange(time=T_EG, rate=m_AF_EU, matrix_index=(1, 0)),
+            msprime.PopulationParametersChange(
+                time=T_EG, growth_rate=r_EU0, initial_size=N_EU1, population_id=1
+            ),
+            msprime.PopulationParametersChange(
+                time=T_EG, growth_rate=0, initial_size=N_AF1, population_id=0
+            ),
+            msprime.MigrationRateChange(time=T_EU0, rate=m_AF_B, matrix_index=(0, 1)),
+            msprime.MigrationRateChange(time=T_EU0, rate=m_AF_B, matrix_index=(1, 0)),
+            msprime.PopulationParametersChange(
+                time=T_EU0, initial_size=N_B, growth_rate=0, population_id=1
+            ),
+            msprime.MassMigration(time=T_OOA, source=1, destination=0, proportion=1.0),
+            msprime.MigrationRateChange(time=T_OOA, rate=0),
+            msprime.PopulationParametersChange(
+                time=T_AF, initial_size=N_A, population_id=0
+            ),
+        ],
+    )
+
+
+_species.add_demographic_model(_ooa_2())
+
+
+def _african():
+    id = "Africa_1T12"
+    description = "African population"
+    long_description = """
+        The model is a simplification of the two population Tennesen et al.
+        model with the European-American population removed so that we are
+        modeling the African population in isolation.
+    """
+    populations = [
+        stdpopsim.Population(id="AFR", description="African"),
+    ]
+    citations = [_tennessen_et_al]
+
+    generation_time = 25
+
+    T_AF = 148e3 / generation_time
+    T_EG = 5115 / generation_time
+
+    # Growth rate
+    r_AF = 0.0166
+
+    # population sizes
+    N_A = 7310
+    N_AF1 = 14474
+
+    # present Ne
+    N_AF = N_AF1 / math.exp(-r_AF * T_EG)
+
+    return stdpopsim.DemographicModel(
+        id=id,
+        description=description,
+        long_description=long_description,
+        populations=populations,
+        citations=citations,
+        generation_time=generation_time,
+        population_configurations=[
+            msprime.PopulationConfiguration(
+                initial_size=N_AF, growth_rate=r_AF, metadata=populations[0].asdict()
+            ),
+        ],
+        demographic_events=[
+            msprime.PopulationParametersChange(
+                time=T_EG, growth_rate=0, initial_size=N_AF1, population_id=0
+            ),
+            msprime.PopulationParametersChange(
+                time=T_AF, initial_size=N_A, population_id=0
+            ),
+        ],
+    )
+
+
+_species.add_demographic_model(_african())
+
+
+def _america():
+    id = "AmericanAdmixture_4B11"
+    description = "American admixture"
+    long_description = """
+        Demographic model for American admixture, taken from Browning et al. 2011.
+        This model extends the Gravel et al. (2011) model of African/European/Asian
+        demographic history to simulate an admixed population with admixture
+        occurring 12 generations ago. The admixed population had an initial size
+        of 30,000 and grew at a rate of 5% per generation, with 1/6 of the
+        population of African ancestry, 1/3 European, and 1/2 Asian.
+    """
+    populations = [
+        stdpopsim.Population(id="AFR", description="Contemporary African population"),
+        stdpopsim.Population(id="EUR", description="Contemporary European population"),
+        stdpopsim.Population(id="ASIA", description="Contemporary Asian population"),
+        stdpopsim.Population(id="ADMIX", description="Modern admixed population"),
+    ]
+
+    citations = [
+        stdpopsim.Citation(
+            author="Browning et al.",
+            year=2011,
+            doi="http://dx.doi.org/10.1371/journal.pgen.1007385",
+            reasons={stdpopsim.CiteReason.DEM_MODEL},
+        )
+    ]
+
+    generation_time = 25
+
+    # Model code was ported from Supplementary File 1.
+    N0 = 7310  # initial population size
+    Thum = 5920  # time (gens) of advent of modern humans
+    Naf = 14474  # size of african population
+    Tooa = 2040  # number of generations back to Out of Africa
+    Nb = 1861  # size of out of Africa population
+    mafb = 1.5e-4  # migration rate Africa and Out-of-Africa
+    Teu = 920  # number generations back to Asia-Europe split
+    Neu = 1032  # bottleneck population sizes
+    Nas = 554
+    mafeu = 2.5e-5  # mig. rates
+    mafas = 7.8e-6
+    meuas = 3.11e-5
+    reu = 0.0038  # growth rate per generation in Europe
+    ras = 0.0048  # growth rate per generation in Asia
+    Tadmix = 12  # time of admixture
+    Nadmix = 30000  # initial size of admixed population
+    radmix = 0.05  # growth rate of admixed population
+    # pop0 is Africa, pop1 is Europe, pop2 is Asia,  pop3 is admixed
+
+    population_configurations = [
+        msprime.PopulationConfiguration(
+            initial_size=Naf, growth_rate=0.0, metadata=populations[0].asdict()
+        ),
+        msprime.PopulationConfiguration(
+            initial_size=Neu * math.exp(reu * Teu),
+            growth_rate=reu,
+            metadata=populations[1].asdict(),
+        ),
+        msprime.PopulationConfiguration(
+            initial_size=Nas * math.exp(ras * Teu),
+            growth_rate=ras,
+            metadata=populations[2].asdict(),
+        ),
+        msprime.PopulationConfiguration(
+            initial_size=Nadmix * math.exp(radmix * Tadmix),
+            growth_rate=radmix,
+            metadata=populations[3].asdict(),
+        ),
+    ]
+
+    migration_matrix = [
+        [0, mafeu, mafas, 0],
+        [mafeu, 0, meuas, 0],
+        [mafas, meuas, 0, 0],
+        [0, 0, 0, 0],
+    ]
+    # Admixture event, 1/6 Africa, 2/6 Europe, 3/6 Asia
+    admixture_event = [
+        msprime.MassMigration(
+            time=Tadmix, source=3, destination=0, proportion=1.0 / 6.0
+        ),
+        msprime.MassMigration(
+            time=Tadmix, source=3, destination=1, proportion=2.0 / 5.0
+        ),
+        msprime.MassMigration(time=Tadmix, source=3, destination=2, proportion=1.0),
+    ]
+    # Asia and Europe split
+    eu_event = [
+        msprime.MigrationRateChange(time=Teu, rate=0.0),
+        msprime.MassMigration(time=Teu, source=2, destination=1, proportion=1.0),
+        msprime.PopulationParametersChange(
+            time=Teu, initial_size=Nb, growth_rate=0.0, population_id=1
+        ),
+        msprime.MigrationRateChange(time=Teu, rate=mafb, matrix_index=(0, 1)),
+        msprime.MigrationRateChange(time=Teu, rate=mafb, matrix_index=(1, 0)),
+    ]
+    # Out of Africa event
+    ooa_event = [
+        msprime.MigrationRateChange(time=Tooa, rate=0.0),
+        msprime.MassMigration(time=Tooa, source=1, destination=0, proportion=1.0),
+    ]
+    # initial population size
+    init_event = [
+        msprime.PopulationParametersChange(time=Thum, initial_size=N0, population_id=0)
+    ]
+    demographic_events = admixture_event + eu_event + ooa_event + init_event
+
+    return stdpopsim.DemographicModel(
+        id=id,
+        description=description,
+        long_description=long_description,
+        populations=populations,
+        citations=citations,
+        generation_time=generation_time,
+        population_configurations=population_configurations,
+        migration_matrix=migration_matrix,
+        demographic_events=demographic_events,
+    )
+
+
+_species.add_demographic_model(_america())
+
+
+def _ooa_archaic():
+    id = "OutOfAfricaArchaicAdmixture_5R19"
+    description = "Three population out-of-Africa with archaic admixture"
+    long_description = """
+        The three population out-of-African model popularized by Gutenkunst et al. (2009)
+        and augmented by archaic contributions to both Eurasian and African populations.
+        Two archaic populations split early in human history, before the African
+        expansion, and contribute to Eurasian populations (putative Neanderthal branch)
+        and to the African branch (a deep diverging branch within Africa). Admixture
+        is modeled as symmetric migration between the archaic and modern human branches,
+        with contribution ending at a given time in the past.
+    """
+    populations = [
+        _yri_population,
+        _ceu_population,
+        _chb_population,
+        stdpopsim.Population(
+            "Neanderthal", "Putative Neanderthals", sampling_time=None
+        ),
+        stdpopsim.Population(
+            "ArchaicAFR", "Putative Archaic Africans", sampling_time=None
+        ),
+    ]
+    citations = [
+        stdpopsim.Citation(
+            author="Ragsdale and Gravel",
+            year=2019,
+            doi="https://doi.org/10.1371/journal.pgen.1008204",
+            reasons={stdpopsim.CiteReason.DEM_MODEL},
+        )
+    ]
+
+    # First we set out the maximum likelihood values of the various parameters
+    # given in Table 1 (under archaic admixture).
+    N_0 = 3600
+    N_YRI = 13900
+    N_B = 880
+    N_CEU0 = 2300
+    N_CHB0 = 650
+
+    # Times are provided in years, so we convert into generations.
+    # In the published model, the authors used a generation time of 29 years to
+    # convert from genetic to physical units
+    generation_time = 29
+
+    T_AF = 300e3 / generation_time
+    T_B = 60.7e3 / generation_time
+    T_EU_AS = 36.0e3 / generation_time
+    T_arch_afr_split = 499e3 / generation_time
+    T_arch_afr_mig = 125e3 / generation_time
+    T_nean_split = 559e3 / generation_time
+    T_arch_adm_end = 18.7e3 / generation_time
+
+    # We need to work out the starting (diploid) population sizes based on
+    # the growth rates provided for these two populations
+    r_CEU = 0.00125
+    r_CHB = 0.00372
+    N_CEU = N_CEU0 / math.exp(-r_CEU * T_EU_AS)
+    N_CHB = N_CHB0 / math.exp(-r_CHB * T_EU_AS)
+
+    # Migration rates during the various epochs.
+    m_AF_B = 52.2e-5
+    m_YRI_CEU = 2.48e-5
+    m_YRI_CHB = 0e-5
+    m_CEU_CHB = 11.3e-5
+    m_AF_arch_af = 1.98e-5
+    m_OOA_nean = 0.825e-5
+
+    # Population IDs correspond to their indexes in the population
+    # configuration array. Therefore, we have 0=YRI, 1=CEU and 2=CHB
+    # initially.
+    # We also have two archaic populations, putative Neanderthals and
+    # archaicAfrican, which are population indices 3=Nean and 4=arch_afr.
+    # Their sizes are equal to the ancestral reference population size N_0.
+    population_configurations = [
+        msprime.PopulationConfiguration(
+            initial_size=N_YRI, metadata=populations[0].asdict()
+        ),
+        msprime.PopulationConfiguration(
+            initial_size=N_CEU, growth_rate=r_CEU, metadata=populations[1].asdict()
+        ),
+        msprime.PopulationConfiguration(
+            initial_size=N_CHB, growth_rate=r_CHB, metadata=populations[2].asdict()
+        ),
+        msprime.PopulationConfiguration(
+            initial_size=N_0, metadata=populations[3].asdict()
+        ),
+        msprime.PopulationConfiguration(
+            initial_size=N_0, metadata=populations[4].asdict()
+        ),
+    ]
+    migration_matrix = [  # noqa
+        [0, m_YRI_CEU, m_YRI_CHB, 0, 0],  # noqa
+        [m_YRI_CEU, 0, m_CEU_CHB, 0, 0],  # noqa
+        [m_YRI_CHB, m_CEU_CHB, 0, 0, 0],  # noqa
+        [0, 0, 0, 0, 0],  # noqa
+        [0, 0, 0, 0, 0],  # noqa
+    ]  # noqa
+    demographic_events = [
+        # first event is migration turned on between modern and archaic humans
+        msprime.MigrationRateChange(
+            time=T_arch_adm_end, rate=m_AF_arch_af, matrix_index=(0, 4)
+        ),
+        msprime.MigrationRateChange(
+            time=T_arch_adm_end, rate=m_AF_arch_af, matrix_index=(4, 0)
+        ),
+        msprime.MigrationRateChange(
+            time=T_arch_adm_end, rate=m_OOA_nean, matrix_index=(1, 3)
+        ),
+        msprime.MigrationRateChange(
+            time=T_arch_adm_end, rate=m_OOA_nean, matrix_index=(3, 1)
+        ),
+        msprime.MigrationRateChange(
+            time=T_arch_adm_end, rate=m_OOA_nean, matrix_index=(2, 3)
+        ),
+        msprime.MigrationRateChange(
+            time=T_arch_adm_end, rate=m_OOA_nean, matrix_index=(3, 2)
+        ),
+        # CEU and CHB merge into B with rate changes at T_EU_AS
+        msprime.MassMigration(time=T_EU_AS, source=2, destination=1, proportion=1.0),
+        msprime.MigrationRateChange(time=T_EU_AS, rate=0),
+        msprime.MigrationRateChange(time=T_EU_AS, rate=m_AF_B, matrix_index=(0, 1)),
+        msprime.MigrationRateChange(time=T_EU_AS, rate=m_AF_B, matrix_index=(1, 0)),
+        msprime.MigrationRateChange(
+            time=T_EU_AS, rate=m_AF_arch_af, matrix_index=(0, 4)
+        ),
+        msprime.MigrationRateChange(
+            time=T_EU_AS, rate=m_AF_arch_af, matrix_index=(4, 0)
+        ),
+        msprime.MigrationRateChange(time=T_EU_AS, rate=m_OOA_nean, matrix_index=(1, 3)),
+        msprime.MigrationRateChange(time=T_EU_AS, rate=m_OOA_nean, matrix_index=(3, 1)),
+        msprime.PopulationParametersChange(
+            time=T_EU_AS, initial_size=N_B, growth_rate=0, population_id=1
+        ),
+        # Population B merges into YRI at T_B
+        msprime.MassMigration(time=T_B, source=1, destination=0, proportion=1.0),
+        msprime.MigrationRateChange(time=T_B, rate=0),
+        msprime.MigrationRateChange(time=T_B, rate=m_AF_arch_af, matrix_index=(0, 4)),
+        msprime.MigrationRateChange(time=T_B, rate=m_AF_arch_af, matrix_index=(4, 0)),
+        # Beginning of migration between African and archaic African populations
+        msprime.MigrationRateChange(time=T_arch_afr_mig, rate=0),
+        # Size changes to N_0 at T_AF
+        msprime.PopulationParametersChange(
+            time=T_AF, initial_size=N_0, population_id=0
+        ),
+        # Archaic African merges with moderns
+        msprime.MassMigration(
+            time=T_arch_afr_split, source=4, destination=0, proportion=1.0
+        ),
+        # Neanderthal merges with moderns
+        msprime.MassMigration(
+            time=T_nean_split, source=3, destination=0, proportion=1.0
+        ),
+    ]
+
+    return stdpopsim.DemographicModel(
+        id=id,
+        description=description,
+        long_description=long_description,
+        populations=populations,
+        citations=citations,
+        generation_time=generation_time,
+        population_configurations=population_configurations,
+        migration_matrix=migration_matrix,
+        demographic_events=demographic_events,
+    )
+
+
+_species.add_demographic_model(_ooa_archaic())
+
+
+def _zigzag():
+    id = "Zigzag_1S14"
+    description = "Periodic growth and decline."
+    long_description = """
+        A validation model used by Schiffels and Durbin (2014) and Terhorst and
+        Terhorst, Kamm, and Song (2017) with periods of exponential growth and
+        decline in a single population.
+        """
+    populations = [
+        stdpopsim.Population("generic", "Generic expanding and contracting population"),
+    ]
+    citations = [
+        stdpopsim.Citation(
+            author="Schiffels and Durbin",
+            year=2014,
+            doi="https://doi.org/10.1038/ng.3015",
+            reasons={stdpopsim.CiteReason.DEM_MODEL},
+        )
+    ]
+
+    generation_time = 30
+    N0 = 5 * 14312
+
+    g_1 = 0.02302578256
+    t_1 = 33.333334976
+
+    g_2 = -0.0057564631
+    t_2 = 133.3334544
+
+    g_3 = 0.00143911578
+    t_3 = 533.33324512
+
+    g_4 = -0.0003597785
+    t_4 = 2133.3352704
+
+    g_5 = 8.9944801565e-5
+    t_5 = 8533.329632
+
+    n_ancient = N0 / 10
+    t_ancient = 34133.318528
+
+    population_configurations = [
+        msprime.PopulationConfiguration(
+            initial_size=N0, metadata=populations[0].asdict()
+        )
+    ]
+
+    demographic_events = [
+        msprime.PopulationParametersChange(time=t_1, growth_rate=g_1, population_id=0),
+        msprime.PopulationParametersChange(time=t_2, growth_rate=g_2, population_id=0),
+        msprime.PopulationParametersChange(time=t_3, growth_rate=g_3, population_id=0),
+        msprime.PopulationParametersChange(time=t_4, growth_rate=g_4, population_id=0),
+        msprime.PopulationParametersChange(time=t_5, growth_rate=g_5, population_id=0),
+        msprime.PopulationParametersChange(
+            time=t_ancient, initial_size=n_ancient, growth_rate=0, population_id=0
+        ),
+    ]
+
+    return stdpopsim.DemographicModel(
+        id=id,
+        description=description,
+        long_description=long_description,
+        populations=populations,
+        citations=citations,
+        generation_time=generation_time,
+        population_configurations=population_configurations,
+        demographic_events=demographic_events,
+    )
+
+
+_species.add_demographic_model(_zigzag())
+
+
+def _kamm_ancient_eurasia():
+    id = "AncientEurasia_9K19"
+    description = "Multi-population model of ancient Eurasia"
+    long_description = """
+        This is the best-fitting model of a history of
+        multiple ancient and present-day human populations
+        sampled across Eurasia over the past 120,000 years.
+        The fitting was performed using momi2 (Kamm et al. 2019),
+        which uses the multi-population site-frequency spectrum
+        as input data. The model includes a ghost admixture event
+        (from unsampled basal Eurasians into early European
+        farmers), and two admixture events where the source is
+        approximately well-known (from Neanderthals into
+        Non-Africans and from Western European hunter-gatherers
+        into modern Sardinians. There are three present-day
+        populations: Sardinians, Han Chinese and African Mbuti.
+        Additionally, there are several ancient samples
+        obtained from fossils dated at different times in
+        the past: the Altai Neanderthal (Prufer et al. 2014),
+        a Mesolithic hunter-gatherer (Lazaridis et al. 2014),
+        a Neolithic early European sample (Lazaridis et al. 2014),
+        and two Palaeolithic modern humans from Siberia - MA1
+        (Raghavan et al. 2014) and  Ust'Ishim (Fu et al. 2014).
+        All the ancient samples are represented by a single diploid
+        genome.
+    """
+    # Sampling times are assuming 25 years per generation
+    populations = [
+        stdpopsim.Population(
+            id="Mbuti", description="Present-day African Mbuti", sampling_time=0
+        ),
+        # LBK: 8,000 years ago
+        stdpopsim.Population(
+            id="LBK", description="Early European farmer (EEF)", sampling_time=320
+        ),
+        stdpopsim.Population(
+            id="Sardinian", description="Present-day Sardinian", sampling_time=0
+        ),
+        # Loschbour: 7,500 years ago
+        stdpopsim.Population(
+            id="Loschbour",
+            description="Western hunter-gatherer (WHG)",
+            sampling_time=300,
+        ),
+        # MA1: 24,000 years ago
+        stdpopsim.Population(
+            id="MA1", description="Upper Palaeolithic MAl'ta culture", sampling_time=960
+        ),
+        stdpopsim.Population(
+            id="Han", description="Present-day Han Chinese", sampling_time=0
+        ),
+        # Ust Ishim: 45,000 years ago
+        stdpopsim.Population(
+            id="UstIshim",
+            description="early Palaeolithic Ust'-Ishim",
+            sampling_time=1800,
+        ),
+        # Altai Neanderthal: 50,000 years ago
+        stdpopsim.Population(
+            id="Neanderthal",
+            description="Altai Neanderthal from Siberia",
+            sampling_time=2000,
+        ),
+        stdpopsim.Population(
+            id="BasalEurasian", description="Basal Eurasians", sampling_time=None
+        ),
+    ]
+    citations = [
+        stdpopsim.Citation(
+            author="Kamm et al.",
+            year=2019,
+            doi="https://doi.org/10.1080/01621459.2019.1635482",
+            reasons={stdpopsim.CiteReason.DEM_MODEL},
+        )
+    ]
+
+    # Times are provided in years, so we convert into generations.
+    generation_time = 25
+    # Mutation_rate in Kamm et al. = 1.22e-8
+    # Effective population sizes
+    N_Losch = 1920
+    N_Mbu = 17300
+    N_Mbu_Losch = 29100
+    N_Han = 6300
+    N_Han_Losch = 2340
+    N_Nean_Losch = 18200
+    N_Nean = 86.9
+    N_LBK = 75.7
+    N_Sard = 15000
+    N_Sard_LBK = 12000
+    # Table A.1 has Altai at 50,000 years ago
+    t_NeaPopSizeChange = 50000 / generation_time
+    # Unknown but suspected parameters...
+    N_Basal = N_Losch
+    N_MA1 = N_Losch
+    N_Ust = N_Losch
+    # Population split times
+    t_Mbu_Losch = 95800 / generation_time
+    t_Han_Losch = 50400 / generation_time
+    t_Ust_Losch = 51500 / generation_time
+    t_Nean_Losch = 696000 / generation_time
+    t_MA1_Losch = 44900 / generation_time
+    t_LBK_Losch = 37700 / generation_time
+    t_Basal_Losch = 79800 / generation_time
+    t_Sard_LBK = 7690 / generation_time
+    # Given that we're using best model estimate,
+    # ghost WHG is directly descended from Loschbour,
+    # so this parameter is not used
+    # t_GhostWHG_Losch = 1560 / generation_time
+    # Admixture times
+    t_Nean_to_Eurasian = 56800 / generation_time
+    t_Basal_to_EEF = 33700 / generation_time
+    t_GhostWHG_to_Sard = 1230 / generation_time
+    t_NeanGrowth = t_Mbu_Losch - t_NeaPopSizeChange
+    logdiffNeanGrowth = math.log(N_Nean / N_Nean_Losch)
+    r_NeanGrowth = logdiffNeanGrowth / t_NeanGrowth
+    p_Nean_to_Eurasian = 0.0296
+    p_Basal_to_EEF = 0.0936
+    p_GhostWHG_to_Sard = 0.0317
+    # Population IDs: Mbuti = 0; LBK = 1;
+    # Sardinian = 2; Loschbour = 3; MA1 = 4;
+    # Han = 5; Ust Ishim = 6; Neanderthal = 7;
+    # Basal Eurasian = 8
+    population_configurations = [
+        msprime.PopulationConfiguration(
+            initial_size=N_Mbu, metadata=populations[0].asdict()
+        ),
+        msprime.PopulationConfiguration(
+            initial_size=N_LBK, metadata=populations[1].asdict()
+        ),
+        msprime.PopulationConfiguration(
+            initial_size=N_Sard, metadata=populations[2].asdict()
+        ),
+        msprime.PopulationConfiguration(
+            initial_size=N_Losch, metadata=populations[3].asdict()
+        ),
+        msprime.PopulationConfiguration(
+            initial_size=N_MA1, metadata=populations[4].asdict()
+        ),
+        msprime.PopulationConfiguration(
+            initial_size=N_Han, metadata=populations[5].asdict()
+        ),
+        msprime.PopulationConfiguration(
+            initial_size=N_Ust, metadata=populations[6].asdict()
+        ),
+        msprime.PopulationConfiguration(
+            initial_size=N_Nean, metadata=populations[7].asdict()
+        ),
+        msprime.PopulationConfiguration(
+            initial_size=N_Basal, metadata=populations[8].asdict()
+        ),
+    ]
+    demographic_events = [
+        # Sardinian receives admixture from Loschbour / WHG
+        msprime.MassMigration(
+            time=t_GhostWHG_to_Sard,
+            source=2,
+            destination=3,
+            proportion=p_GhostWHG_to_Sard,
+        ),
+        # Sardinian merges into LBK / EEF
+        # Now pop 1: Sardinian-LBK ancestral pop
+        msprime.MassMigration(time=t_Sard_LBK, source=2, destination=1, proportion=1.0),
+        # Sardinian-LBK ancestral pop size change
+        msprime.PopulationParametersChange(
+            time=t_Sard_LBK, initial_size=N_Sard_LBK, population_id=1
+        ),
+        # LBK / EEF receives admixture from Basal Eurasians
+        msprime.MassMigration(
+            time=t_Basal_to_EEF, source=1, destination=8, proportion=p_Basal_to_EEF
+        ),
+        # LBK / EEF merges into Loschbour
+        msprime.MassMigration(
+            time=t_LBK_Losch, source=1, destination=3, proportion=1.0
+        ),
+        # MA1 merges into Loschbour
+        msprime.MassMigration(
+            time=t_MA1_Losch, source=4, destination=3, proportion=1.0
+        ),
+        # Neanderthal start change in population size
+        msprime.PopulationParametersChange(
+            time=t_NeaPopSizeChange,
+            initial_size=N_Nean,
+            growth_rate=r_NeanGrowth,
+            population_id=7,
+        ),
+        # Han merges into Loschbour
+        msprime.MassMigration(
+            time=t_Han_Losch, source=5, destination=3, proportion=1.0
+        ),
+        # Change in population size in Han-Losch ancestral pop
+        msprime.PopulationParametersChange(
+            time=t_Han_Losch, initial_size=N_Han_Losch, population_id=3
+        ),
+        # UstIshim merges into Loschbour
+        msprime.MassMigration(
+            time=t_Ust_Losch, source=6, destination=3, proportion=1.0
+        ),
+        # Loschbour / Non-Africans receive admixture from Neanderthals
+        msprime.MassMigration(
+            time=t_Nean_to_Eurasian,
+            source=3,
+            destination=7,
+            proportion=p_Nean_to_Eurasian,
+        ),
+        # Basal Eurasians merge into Loschbour / Non-Africans
+        msprime.MassMigration(
+            time=t_Basal_Losch, source=8, destination=3, proportion=1.0
+        ),
+        # Mbuti merge into Loschbour / Non-Africans
+        msprime.MassMigration(
+            time=t_Mbu_Losch, source=0, destination=3, proportion=1.0
+        ),
+        # Change in population size in Mbuti-Losch ancestral pop
+        msprime.PopulationParametersChange(
+            time=t_Mbu_Losch, initial_size=N_Mbu_Losch, population_id=3
+        ),
+        # Change in population size in Neanderthal, growth rate 0
+        msprime.PopulationParametersChange(
+            time=t_Mbu_Losch, initial_size=N_Nean_Losch, growth_rate=0, population_id=7
+        ),
+        # Neanderthal merge into Loschbour / modern humans
+        msprime.MassMigration(
+            time=t_Nean_Losch, source=7, destination=3, proportion=1.0
+        ),
+        # Ancestral hominin population size change
+        msprime.PopulationParametersChange(
+            time=t_Nean_Losch, initial_size=N_Nean_Losch, population_id=3
+        ),
+    ]
+
+    return stdpopsim.DemographicModel(
+        id=id,
+        description=description,
+        long_description=long_description,
+        populations=populations,
+        citations=citations,
+        generation_time=generation_time,
+        population_configurations=population_configurations,
+        demographic_events=demographic_events,
+    )
+
+
+_species.add_demographic_model(_kamm_ancient_eurasia())
+
+
+def _papuans_10j19():
+    id = "PapuansOutOfAfrica_10J19"
+    description = "Out-of-Africa with archaic admixture into Papuans"
+    long_description = """
+        A ten population model of out-of-Africa, including two pulses of
+        Denisovan admixture into Papuans, and several pulses of Neandertal
+        admixture into non-Africans.
+        Most parameters are from Jacobs et al. (2019), Table S5 and Figure S5.
+        This model is an extension of one from Malaspinas et al. (2016), thus
+        some parameters are inherited from there.
+        """
+
+    # sampling times
+    T_DenA = 2058
+    T_NeaA = 2612
+
+    populations = [
+        # humans
+        _yri_population,
+        _ceu_population,
+        _chb_population,
+        stdpopsim.Population("Papuan", "Papuans from Indonesia and New Guinea"),
+        # archaics
+        stdpopsim.Population(
+            "DenA", "Altai Denisovan (sampling) lineage", sampling_time=T_DenA
+        ),
+        stdpopsim.Population(
+            "NeaA", "Altai Neandertal (sampling) lineage", sampling_time=T_NeaA
+        ),
+        stdpopsim.Population(
+            "Den1", "Denisovan D1 (introgressing) lineage", sampling_time=None
+        ),
+        stdpopsim.Population(
+            "Den2", "Denisovan D2 (introgressing) lineage", sampling_time=None
+        ),
+        stdpopsim.Population(
+            "Nea1", "Neandertal N1 (introgressing) lineage", sampling_time=None
+        ),
+        stdpopsim.Population("Ghost", "Out-of-Africa lineage", sampling_time=None),
+    ]
+    pop = {p.id: i for i, p in enumerate(populations)}
+
+    citations = [
+        stdpopsim.Citation(
+            author="Jacobs et al.",
+            year=2019,
+            doi="https://doi.org/10.1016/j.cell.2019.02.035",
+            reasons={stdpopsim.CiteReason.DEM_MODEL},
+        ),
+        stdpopsim.Citation(
+            author="Malaspinas et al.",
+            year=2016,
+            doi="https://doi.org/10.1038/nature18299",
+            reasons={stdpopsim.CiteReason.DEM_MODEL},
+        ),
+    ]
+
+    # Inherited from Malaspinas et al., which gives the following refs:
+    generation_time = 29  # Fenner 2005
+    # mutation_rate = 1.25e-8  # per gen per site, Scally & Durbin 2012
+
+    N_YRI = 48433
+    N_Ghost = 8516
+    N_CEU = 6962
+    N_CHB = 9025
+    N_Papuan = 8834
+
+    N_DenA = 5083
+    N_Den1 = 13249
+    N_Den2 = 13249
+    N_NeaA = 826
+    N_Nea1 = 13249
+
+    pop_meta = {p.id: p.asdict() for p in populations}
+    population_configurations = [
+        msprime.PopulationConfiguration(initial_size=N_YRI, metadata=pop_meta["YRI"]),
+        msprime.PopulationConfiguration(initial_size=N_CEU, metadata=pop_meta["CEU"]),
+        msprime.PopulationConfiguration(initial_size=N_CHB, metadata=pop_meta["CHB"]),
+        msprime.PopulationConfiguration(
+            initial_size=N_Papuan, metadata=pop_meta["Papuan"]
+        ),
+        msprime.PopulationConfiguration(initial_size=N_DenA, metadata=pop_meta["DenA"]),
+        msprime.PopulationConfiguration(initial_size=N_NeaA, metadata=pop_meta["NeaA"]),
+        msprime.PopulationConfiguration(initial_size=N_Den1, metadata=pop_meta["Den1"]),
+        msprime.PopulationConfiguration(initial_size=N_Den2, metadata=pop_meta["Den2"]),
+        msprime.PopulationConfiguration(initial_size=N_Nea1, metadata=pop_meta["Nea1"]),
+        msprime.PopulationConfiguration(
+            initial_size=N_Ghost, metadata=pop_meta["Ghost"]
+        ),
+    ]
+    assert len(populations) == len(population_configurations)
+
+    # initial migrations
+    m_Ghost_Afr = 1.79e-4
+    m_Ghost_EU = 4.42e-4
+    m_EU_AS = 3.14e-5
+    m_AS_Papuan = 5.72e-5
+    # older migrations
+    m_Eurasian_Papuan = 5.72e-4
+    m_Eurasian_Ghost = 4.42e-4
+
+    migration_matrix = [[0] * len(populations) for _ in range(len(populations))]
+    migration_matrix[pop["Ghost"]][pop["YRI"]] = m_Ghost_Afr
+    migration_matrix[pop["YRI"]][pop["Ghost"]] = m_Ghost_Afr
+    migration_matrix[pop["Ghost"]][pop["CEU"]] = m_Ghost_EU
+    migration_matrix[pop["CEU"]][pop["Ghost"]] = m_Ghost_EU
+    migration_matrix[pop["CEU"]][pop["CHB"]] = m_EU_AS
+    migration_matrix[pop["CHB"]][pop["CEU"]] = m_EU_AS
+    migration_matrix[pop["CHB"]][pop["Papuan"]] = m_AS_Papuan
+    migration_matrix[pop["Papuan"]][pop["CHB"]] = m_AS_Papuan
+
+    # splits
+    T_EU_AS_split = 1293
+    T_Eurasian_Ghost_split = 1758
+    T_Papuan_Ghost_split = 1784
+    T_Ghost_Afr_split = 2218
+    T_NeaA_Nea1_split = 3375
+    T_DenA_Den1_split = 9750
+    T_DenA_Den2_split = 12500
+    T_Den_Nea_split = 15090
+    T_Afr_Archaic_split = 20225
+
+    # bottlenecks
+    Tb_Eurasia = 1659
+    Tb_Papua = 1685
+    Tb_Ghost = 2119
+    Nb_Eurasia = 2231
+    Nb_Papua = 243
+    Nb_Ghost = 1394
+
+    # internal branches
+    N_EU_AS = 12971
+    N_Ghost_Afr = 41563
+    N_NeaA_Nea1 = 13249
+    N_Den_Anc = 100  # S10.i p. 31/45
+    N_DenA_Den1 = N_Den_Anc
+    N_DenA_Den2 = N_Den_Anc
+    N_Den_Nea = 13249
+    N_Afr_Archaic = 32671
+
+    # admixture pulses
+    m_Den_Papuan = 0.04
+    p1 = 0.55  # S10.i p. 31/45
+    m_Den1_Papuan = p1 * m_Den_Papuan
+    m_Den2_Papuan = (1 - p1) * m_Den_Papuan
+    m_Nea1_Ghost = 0.024
+    m_Nea1_Eurasian = 0.011
+    m_Nea1_Papuan = 0.002
+    m_Nea1_AS = 0.002
+
+    T_Nea1_Ghost_mig = 1853
+    T_Nea1_Eurasian_mig = 1566
+    T_Nea1_Papuan_mig = 1412
+    T_Nea1_AS_mig = 883
+    # Fig. 4B, and S10.h p. 30/45
+    T_Den1_Papuan_mig = 29.8e3 / generation_time
+    T_Den2_Papuan_mig = 45.7e3 / generation_time
+
+    demographic_events = [
+        # human lineage splits
+        msprime.MassMigration(
+            time=T_EU_AS_split, source=pop["CEU"], destination=pop["CHB"]
+        ),
+        msprime.PopulationParametersChange(
+            time=T_EU_AS_split, initial_size=N_EU_AS, population_id=pop["CHB"]
+        ),
+        msprime.MassMigration(
+            time=T_Eurasian_Ghost_split, source=pop["CHB"], destination=pop["Ghost"]
+        ),
+        msprime.MassMigration(
+            time=T_Papuan_Ghost_split, source=pop["Papuan"], destination=pop["Ghost"]
+        ),
+        msprime.MassMigration(
+            time=T_Ghost_Afr_split, source=pop["Ghost"], destination=pop["YRI"]
+        ),
+        msprime.PopulationParametersChange(
+            time=T_Ghost_Afr_split, initial_size=N_Ghost_Afr, population_id=pop["YRI"]
+        ),
+        # archaic lineage splits
+        msprime.MassMigration(
+            time=T_DenA_Den1_split, source=pop["Den1"], destination=pop["DenA"]
+        ),
+        msprime.PopulationParametersChange(
+            time=T_DenA_Den1_split, initial_size=N_DenA_Den1, population_id=pop["DenA"]
+        ),
+        msprime.MassMigration(
+            time=T_DenA_Den2_split, source=pop["Den2"], destination=pop["DenA"]
+        ),
+        msprime.PopulationParametersChange(
+            time=T_DenA_Den2_split, initial_size=N_DenA_Den2, population_id=pop["DenA"]
+        ),
+        msprime.MassMigration(
+            time=T_NeaA_Nea1_split, source=pop["Nea1"], destination=pop["NeaA"]
+        ),
+        msprime.PopulationParametersChange(
+            time=T_NeaA_Nea1_split, initial_size=N_NeaA_Nea1, population_id=pop["NeaA"]
+        ),
+        msprime.MassMigration(
+            time=T_Den_Nea_split, source=pop["NeaA"], destination=pop["DenA"]
+        ),
+        msprime.PopulationParametersChange(
+            time=T_Den_Nea_split, initial_size=N_Den_Nea, population_id=pop["DenA"]
+        ),
+        msprime.MassMigration(
+            time=T_Afr_Archaic_split, source=pop["DenA"], destination=pop["YRI"]
+        ),
+        msprime.PopulationParametersChange(
+            time=T_Afr_Archaic_split,
+            initial_size=N_Afr_Archaic,
+            population_id=pop["YRI"],
+        ),
+        # bottlenecks
+        msprime.PopulationParametersChange(
+            time=Tb_Eurasia, initial_size=Nb_Eurasia, population_id=pop["CHB"]
+        ),
+        msprime.PopulationParametersChange(
+            time=Tb_Papua, initial_size=Nb_Papua, population_id=pop["Papuan"]
+        ),
+        msprime.PopulationParametersChange(
+            time=Tb_Ghost, initial_size=Nb_Ghost, population_id=pop["Ghost"]
+        ),
+        # migration changes
+        msprime.MigrationRateChange(
+            time=T_EU_AS_split, rate=0, matrix_index=(pop["CHB"], pop["CEU"])
+        ),
+        msprime.MigrationRateChange(
+            time=T_EU_AS_split, rate=0, matrix_index=(pop["CEU"], pop["CHB"])
+        ),
+        msprime.MigrationRateChange(
+            time=T_EU_AS_split, rate=0, matrix_index=(pop["Papuan"], pop["CHB"])
+        ),
+        msprime.MigrationRateChange(
+            time=T_EU_AS_split, rate=0, matrix_index=(pop["CHB"], pop["Papuan"])
+        ),
+        msprime.MigrationRateChange(
+            time=T_EU_AS_split, rate=0, matrix_index=(pop["Ghost"], pop["CEU"])
+        ),
+        msprime.MigrationRateChange(
+            time=T_EU_AS_split, rate=0, matrix_index=(pop["CEU"], pop["Ghost"])
+        ),
+        msprime.MigrationRateChange(
+            time=T_EU_AS_split,
+            rate=m_Eurasian_Papuan,
+            matrix_index=(pop["CHB"], pop["Papuan"]),
+        ),
+        msprime.MigrationRateChange(
+            time=T_EU_AS_split,
+            rate=m_Eurasian_Papuan,
+            matrix_index=(pop["Papuan"], pop["CHB"]),
+        ),
+        msprime.MigrationRateChange(
+            time=T_EU_AS_split,
+            rate=m_Eurasian_Ghost,
+            matrix_index=(pop["CHB"], pop["Ghost"]),
+        ),
+        msprime.MigrationRateChange(
+            time=T_EU_AS_split,
+            rate=m_Eurasian_Ghost,
+            matrix_index=(pop["Ghost"], pop["CHB"]),
+        ),
+        # all migrations off
+        msprime.MigrationRateChange(time=Tb_Eurasia, rate=0),
+        # admixture pulses
+        msprime.MassMigration(
+            time=T_Den1_Papuan_mig,
+            proportion=m_Den1_Papuan,
+            source=pop["Papuan"],
+            destination=pop["Den1"],
+        ),
+        msprime.MassMigration(
+            time=T_Den2_Papuan_mig,
+            proportion=m_Den2_Papuan,
+            source=pop["Papuan"],
+            destination=pop["Den2"],
+        ),
+        msprime.MassMigration(
+            time=T_Nea1_Ghost_mig,
+            proportion=m_Nea1_Ghost,
+            source=pop["Ghost"],
+            destination=pop["Nea1"],
+        ),
+        msprime.MassMigration(
+            time=T_Nea1_Eurasian_mig,
+            proportion=m_Nea1_Eurasian,
+            source=pop["CHB"],
+            destination=pop["Nea1"],
+        ),
+        msprime.MassMigration(
+            time=T_Nea1_Papuan_mig,
+            proportion=m_Nea1_Papuan,
+            source=pop["Papuan"],
+            destination=pop["Nea1"],
+        ),
+        msprime.MassMigration(
+            time=T_Nea1_AS_mig,
+            proportion=m_Nea1_AS,
+            source=pop["CHB"],
+            destination=pop["Nea1"],
+        ),
+    ]
+
+    demographic_events.sort(key=lambda x: x.time)
+
+    return stdpopsim.DemographicModel(
+        id=id,
+        description=description,
+        long_description=long_description,
+        populations=populations,
+        citations=citations,
+        generation_time=generation_time,
+        population_configurations=population_configurations,
+        migration_matrix=migration_matrix,
+        demographic_events=demographic_events,
+    )
+
+
+_species.add_demographic_model(_papuans_10j19())
+
+
+def _AJ():
+    id = "AshkSub_7G19"
+    description = "Ashkenazi Jewish with substructure and European admixture"
+    long_description = """
+    This was the best fit model of Ashkenazi Jewish demographic history from
+    Gladstein and Hammer 2019, shown in Figure 1, labeled "Substructure Model".
+    Model choice and parameter estimation were performed with Approximate
+    Bayesian Computation. Parameter values are based on the mode from ABC found
+    in Table S3 of Gladstein and Hammer 2019. In this model, the ancestors of
+    Europeans and Middle Eastern populations diverge. Non-Ashkenazi Jewish
+    populations then diverge from the Middle Eastern population. The Ashkenazi
+    Jews then diverge from the other Jewish populations and experience a
+    substantial reduction in population size and a single pulse of gene flow
+    from Europeans (corresponding to their arrival in Europe). After the gene
+    flow from Europeans to the Ashkenazi Jews, the Ashkenazi Jews split into
+    two groups, the Western and Eastern. Finally, the Western Ashkenazi Jews
+    experience moderate instantaneous population size increase, and the
+    Eastern experience a massive population size increase. In addition to the
+    demographic model Gladstein and Hammer 2019 also incorporated an SNP array
+    ascertainment scheme into the simulation. This demographic model does not
+    include the SNP array ascertainment scheme. It should be noted that
+    Gladstein and Hammer 2019 simulated with a mutation rate of 2.5e-8.
+    """
+    populations = [
+        _yri_population,
+        _chb_population,
+        _ceu_population,
+        stdpopsim.Population(id="ME", description="Middle Eastern"),
+        stdpopsim.Population(id="J", description="non-Ashkenazi Jewish"),
+        stdpopsim.Population(id="WAJ", description="Western Ashkenazi Jewish"),
+        stdpopsim.Population(id="EAJ", description="Eastern Ashkenazi Jewish"),
+    ]
+    citations = [
+        stdpopsim.Citation(
+            author="Gladstein and Hammer",
+            year=2019,
+            doi="https://doi.org/10.1093/molbev/msz047",
+            reasons={stdpopsim.CiteReason.DEM_MODEL},
+        )
+    ]
+
+    generation_time = 25
+
+    # parameter value definitions based on mode from ABC
+    # found in Table S3 of Gladstein and Hammer 2019
+
+    # effective population sizes
+    NANC = 7300  # not inferred. Value taken from Gutenkunst et al. 2009
+    NYRI = 10 ** 4.26
+    NCHB = 10 ** 3.61
+    NCEU = 10 ** 4.52
+    NM = 10 ** 5.64
+    NJ = 10 ** 5.55
+    NAg = 10 ** 3.04
+    NWA = 10 ** 3.82
+    NEA = 10 ** 6.29
+
+    # admixture proportion from European to Ashkenazi Jews.
+    m = 0.17
+
+    # Times in generations
+    Tgrowth = 8800  # not inferred. Value taken from Gutenkunst et al. 2009
+    TAF = 2105
+    Teu_as = 850
+    TEM = 481
+    TMJ = 211
+    TA = 29
+    Tm = 28  # not inferred. For simplicity put at one generation before TA
+    TAEW = 14
+    TAg = 13  # not inferred. For simplicity put at one generation before TAEW
+
+    YRI, CHB, CEU, M, J, WA, EA = 0, 1, 2, 3, 4, 5, 6
+
+    return stdpopsim.DemographicModel(
+        id=id,
+        description=description,
+        long_description=long_description,
+        populations=populations,
+        citations=citations,
+        generation_time=generation_time,
+        population_configurations=[
+            msprime.PopulationConfiguration(
+                initial_size=NYRI, metadata=populations[0].asdict()
+            ),
+            msprime.PopulationConfiguration(
+                initial_size=NCHB, metadata=populations[1].asdict()
+            ),
+            msprime.PopulationConfiguration(
+                initial_size=NCEU, metadata=populations[2].asdict()
+            ),
+            msprime.PopulationConfiguration(
+                initial_size=NM, metadata=populations[3].asdict()
+            ),
+            msprime.PopulationConfiguration(
+                initial_size=NJ, metadata=populations[4].asdict()
+            ),
+            msprime.PopulationConfiguration(
+                initial_size=NWA, metadata=populations[5].asdict()
+            ),
+            msprime.PopulationConfiguration(
+                initial_size=NEA, metadata=populations[6].asdict()
+            ),
+        ],
+        demographic_events=[
+            # instantaneous growth in EA and WA at the same time
+            msprime.PopulationParametersChange(
+                time=TAg, initial_size=NAg, population_id=WA
+            ),
+            msprime.PopulationParametersChange(
+                time=TAg, initial_size=NAg, population_id=EA
+            ),
+            # EA splits from WA
+            msprime.MassMigration(time=TAEW, source=EA, destination=WA, proportion=1.0),
+            # E geneflow into WA ancestor (forward in time)
+            msprime.MassMigration(time=Tm, source=WA, destination=CEU, proportion=m),
+            # WA ancestor splits from J
+            msprime.MassMigration(time=TA, source=WA, destination=J, proportion=1.0),
+            # J splits from M split
+            msprime.MassMigration(time=TMJ, source=J, destination=M, proportion=1.0),
+            # M splits from CEU
+            msprime.MassMigration(time=TEM, source=M, destination=CEU, proportion=1.0),
+            # CEU splits from CHB
+            msprime.MassMigration(
+                time=Teu_as, source=CEU, destination=CHB, proportion=1.0
+            ),
+            # CHB splits from YRI
+            msprime.MassMigration(
+                time=TAF, source=CHB, destination=YRI, proportion=1.0
+            ),
+            # Instantaneous growth in YRI
+            msprime.PopulationParametersChange(
+                time=Tgrowth, initial_size=NANC, population_id=YRI
+            ),
+        ],
+    )
+
+
+_species.add_demographic_model(_AJ())
+
+
+def _ooa_4pop():
+    id = "OutOfAfrica_4J17"
+    description = "4 population out of Africa"
+    long_description = """
+        Demographic model for a four population out-of-Africa history,
+        taken from Jouganous et al. (2017). Parameter values were taken
+        from table 4 in the main text. This model was fit based on joint
+        allele frequecy spectrum (AFS) data from 1000 Genomes exomes
+        from the YRI, CEU, CHB, and JPT poulation samples. The demography
+        follows the previous three-populations out-of-Africa models
+        with an additional population split in Asia leading to the
+        Japanese (JPT) population. Parameter values were estimated with
+        the program Moments assuming a mutation rate of 1.44e-8 and a
+        generation time of 29 years.
+    """
+    populations = [
+        _yri_population,
+        _ceu_population,
+        _chb_population,
+        stdpopsim.Population(
+            id="JPT",
+            description="1000 Genomes JPT (Japanese in Tokyo, Japan)",
+            sampling_time=0,
+        ),
+    ]
+    citations = [
+        stdpopsim.Citation(
+            author="Jouganous et al.",
+            year="2017",
+            doi="https://doi.org/10.1534/genetics.117.200493",
+            reasons={stdpopsim.CiteReason.DEM_MODEL},
+        )
+    ]
+
+    generation_time = 29
+
+    # Parameter values from Table 4 (in bold)
+
+    # Population sizes
+    N_A = 11293  # ancestral YRI population size
+    N_AF = 23721  # YRI population size
+    N_B = 2831  # B population size (out-of-Africa)
+    N_EU0 = 2512  # CEU population size (intial)
+    N_AS0 = 1019  # Asian population size (inital)
+    N_JP0 = 4384  # Japanese population size (inital)
+
+    # growth rates
+    r_EU = 0.0016  # European growth rate
+    r_JP = 0.0129  # Japanese growth rate
+    r_AS = 0.0026  # Asian growth rate
+
+    # migration rates
+    m_AF_B = 16.8e-5  # YRI <--> B (out-of-Africa)
+    m_AF_EU = 1.14e-5  # YRI <--> CEU
+    m_AF_AS = 0.56e-5  # YRI <--> CHB
+    m_EU_AS = 4.75e-5  # European <--> CHB
+    m_CH_JP = 3.3e-5  # CHBHan <--> JPT
+
+    # epochs (specified in years, converted to generations)
+    T_AF = int(357000 / generation_time)  # time ancestral to African
+    T_B = int(119000 / generation_time)  # time of Out of Africa split
+    T_EU_AS = int(46000 / generation_time)  # time of Europe/Asia split
+    T_CH_JP = int(9000 / generation_time)  # time of Han/Japanese split
+
+    # also
+    r_AF = 0  # YRI and B growth rate
+
+    # A and YRI
+    pop0 = msprime.PopulationConfiguration(
+        initial_size=N_AF * math.exp(r_AF * T_AF),
+        growth_rate=r_AF,
+        metadata=populations[0].asdict(),
+    )
+
+    # B and CEU
+    pop1 = msprime.PopulationConfiguration(
+        initial_size=N_EU0 * math.exp(r_EU * T_EU_AS),
+        growth_rate=r_EU,
+        metadata=populations[1].asdict(),
+    )
+
+    # Asian and CHB
+    pop2 = msprime.PopulationConfiguration(
+        initial_size=N_AS0 * math.exp(r_AS * T_EU_AS),
+        growth_rate=r_AS,
+        metadata=populations[2].asdict(),
+    )
+
+    # JPT
+    pop3 = msprime.PopulationConfiguration(
+        initial_size=N_JP0 * math.exp(r_JP * T_CH_JP),
+        growth_rate=r_JP,
+        metadata=populations[3].asdict(),
+    )
+
+    population_configurations = [pop0, pop1, pop2, pop3]
+
+    migration_matrix = [
+        [0, m_AF_EU, m_AF_AS, 0],
+        [m_AF_EU, 0, m_EU_AS, 0],
+        [m_AF_AS, m_EU_AS, 0, m_CH_JP],
+        [0, 0, m_CH_JP, 0],
+    ]
+
+    # Demographic events from most recent -> oldest
+
+    # JPT splits from CHB
+    JPT_event = [
+        # seg migration rate JPT <--> CHB to zero
+        msprime.MigrationRateChange(time=T_CH_JP, rate=0.0, matrix_index=(2, 3)),
+        msprime.MigrationRateChange(time=T_CH_JP, rate=0.0, matrix_index=(3, 2)),
+        # JPT splits from CHB
+        msprime.MassMigration(time=T_CH_JP, source=3, destination=2, proportion=1.0),
+    ]
+
+    # CEU and CHB merge into B
+    AS_event = [
+        # set migration rates to zero
+        msprime.MigrationRateChange(time=T_EU_AS, rate=0.0),
+        # CHB splits from CEU
+        msprime.MassMigration(time=T_EU_AS, source=2, destination=1, proportion=1.0),
+        # change size of B to N_B
+        # set growth rates in B to zero
+        msprime.PopulationParametersChange(
+            time=T_EU_AS, initial_size=N_B, growth_rate=0.0, population_id=1
+        ),
+        # set up migration between YRI and B
+        msprime.MigrationRateChange(time=T_EU_AS, rate=m_AF_B, matrix_index=(0, 1)),
+        msprime.MigrationRateChange(time=T_EU_AS, rate=m_AF_B, matrix_index=(1, 0)),
+    ]
+
+    # Out-of-Africa event
+    EU_event = [
+        # B splits from YRI
+        msprime.MigrationRateChange(time=T_B, rate=0.0),
+        msprime.MassMigration(time=T_B, source=1, destination=0, proportion=1.0),
+    ]
+
+    # Change size of YRI/African population
+    AF_event = [
+        msprime.PopulationParametersChange(time=T_AF, initial_size=N_A, population_id=0)
+    ]
+
+    demographic_events = JPT_event + AS_event + EU_event + AF_event
+
+    return stdpopsim.DemographicModel(
+        id=id,
+        description=description,
+        long_description=long_description,
+        populations=populations,
+        citations=citations,
+        generation_time=generation_time,
+        population_configurations=population_configurations,
+        migration_matrix=migration_matrix,
+        demographic_events=demographic_events,
+    )
+
+
+_species.add_demographic_model(_ooa_4pop())
diff --git a/stdpopsim/catalog/HomSap/genetic_maps.py b/stdpopsim/catalog/HomSap/genetic_maps.py
new file mode 100644
index 000000000..740de6231
--- /dev/null
+++ b/stdpopsim/catalog/HomSap/genetic_maps.py
@@ -0,0 +1,155 @@
+import stdpopsim
+
+
+_hapmap2007 = stdpopsim.Citation(
+    doi="https://doi.org/10.1038/nature06258",
+    year=2007,
+    author="The International HapMap Consortium",
+)
+
+_species = stdpopsim.get_species("HomSap")
+
+_gm = stdpopsim.GeneticMap(
+    species=_species,
+    id="HapMapII_GRCh37",
+    description="HapMap Phase II lifted over to GRCh37",
+    long_description="""
+        This genetic map is from the Phase II Hapmap project
+        and based on 3.1 million genotyped SNPs
+        from 270 individuals across four populations (YRI, CEU, CHB and JPT).
+        Genome wide recombination rates were estimated using LDHat.
+        This version of the HapMap genetic map was lifted over to GRCh37
+        (and adjusted in regions where the genome assembly had rearranged)
+        for use in the 1000 Genomes project. Please see the README file on
+        the 1000 Genomes download site for details of these adjustments.
+        ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/working/20110106_recombination_hotspots
+        """,
+    url=(
+        "https://stdpopsim.s3-us-west-2.amazonaws.com/genetic_maps/"
+        "HomSap/HapmapII_GRCh37_RecombinationHotspots.tar.gz"
+    ),
+    sha256="80f22d9e6cb0e497074ed1bc277e765fa9d8e22f21b2f66c3b10286520f6b68f",
+    file_pattern="genetic_map_GRCh37_chr{id}.txt",
+    citations=[_hapmap2007.because(stdpopsim.CiteReason.GEN_MAP)],
+)
+_species.add_genetic_map(_gm)
+
+_gm = stdpopsim.GeneticMap(
+    species=_species,
+    id="DeCodeSexAveraged_GRCh36",
+    description="Sex averaged map from deCode family study",
+    long_description="""
+        This genetic map is from the deCode study of recombination
+        events in 15,257 parent-offspring pairs from Iceland.
+        289,658 phased autosomal SNPs were used to call recombinations
+        within these families, and recombination rates computed from the
+        density of these events. This is the combined male and female
+        (sex averaged) map. See
+        https://www.decode.com/addendum/ for more details.""",
+    url=(
+        "https://stdpopsim.s3-us-west-2.amazonaws.com/genetic_maps/"
+        "HomSap/decode_2010_sex-averaged_map.tar.gz"
+    ),
+    sha256="a17f491756f863c1cc61240c0ce4f383d8085c68b45edd3fed677bb9692af23d",
+    file_pattern="genetic_map_decode_2010_sex-averaged_chr{id}.txt",
+    citations=[
+        stdpopsim.Citation(
+            year=2010,
+            author="Kong et al",
+            doi="https://doi.org/10.1038/nature09525",
+            reasons={stdpopsim.CiteReason.GEN_MAP},
+        )
+    ],
+)
+_species.add_genetic_map(_gm)
+
+# 26 populations in 1000 Genomes
+for pop, sha256 in zip(
+    [
+        "ACB",
+        "ASW",
+        "BEB",
+        "CDX",
+        "CEU",
+        "CHB",
+        "CHS",
+        "CLM",
+        "ESN",
+        "FIN",
+        "GBR",
+        "GIH",
+        "GWD",
+        "IBS",
+        "ITU",
+        "JPT",
+        "KHV",
+        "LWK",
+        "MSL",
+        "MXL",
+        "PEL",
+        "PJL",
+        "PUR",
+        "STU",
+        "TSI",
+        "YRI",
+    ],
+    """
+588be112b6aac9f82bfb64ea828809051cdf021a627450b2d8c37d7909486903
+e7017a82e92e87f443cd21c7f1c3ce24168fa5efaa48f8448c8e095d13aebcaf
+c6f69e6a847af6f88fbf1c82024af6fd5eeb35c167e3a721d9de6d6ff146a3ba
+37cd4e8ac16775b9c8126dc09ab90840925f47138675d6a5942dd6264d089f6c
+fcfacc1858bff5a64480982bab0821d47c4c56262310dd85290f907b8ea6b081
+3836a7853b2e128b6bdbeb0e5eb1ba582975304a9efa10b742be856efdc520a9
+e0ba174edf944e96f327d799010a5be8ecc110968be6b27fb8f6901d5a413cbb
+969f04541d90b70e68f75c0edcfffeb00bee2cd97cf7ad598439f18a8d3a0885
+6310f3860e4ccfcedddc173df2a03400d654ddd83c78eb88ec92e104d00bf40b
+bb7ac5df072b8b2a4567825c5ee4382b67cae5290f5879a1d7c3f743f244b3f6
+798daacb41ae538017303f103c39cdb88add7d88eb964922fe6c81406d8aaa64
+28e4552e485b55b0bc392485cd389a155fa83ee3556f221b3b3377d2c996cbdf
+f2d57c835ac61763112c15acd9cc8240af7c7fb1d7cfa7543cb65a565da2f62b
+bbc49af3646b47dd5b58b8ae28716dcce695eb4a9cb223e261f45345b1a051b9
+050b965efb0b4c9197a289521406586f4207dd4833be9fd429d118f7d061788a
+e7284218b4bf3e1487b9194c9f243c5dff88f5edbd40247118225a9acc08c76e
+74cd0bc6bc12c2b59f928c6a3d60773b6118625375266ffba0c89da61c376f55
+dd3ef07641ad7c8d3bfa74a98a5820f1b70ea74c3e82103663bbe61931bf8579
+ad5064ee575a10f890c81b9bdbd55128197582db1c9b153b575f5e7fbe542b0d
+f37b50b940816d9c12429a4b71709b18569e77c1ea4c78d5fed19141751822c3
+024ac242da742d006a210cbf77c5b8457735a6690dc78cc01686deafe7cb7beb
+871b518e3b9ca92838e206125e4a1187690018ec5bddeaf39eb2e6f6be5705f1
+00f87ca1ad6430b272251d3ce5fa93b44a0a50d6472237bcb5dbad2dbcd567b8
+fa9e61492de4c160cf987ac6a8cc3a77efdea751f0f5aa9303a294a37f81cbd7
+e3cf06041c6ffceb29fbc617c86966948203e166c3399d036c95915e3d5ebdcd
+2989aecaa086af37cdef9ee2245f12000b3f29c08b5d1a71744578a684fe377f
+    """.strip().splitlines(),
+):
+    _gm = stdpopsim.GeneticMap(
+        species=_species,
+        id="Pyrho{}_GRCh38".format(pop),
+        description="Pyrho population-specific map for {}".format(pop),
+        long_description="""
+        This genetic map was inferred using individuals from the {}
+        population from Phase 3 of the 1000 Genomes Project. Rates
+        were estimated using pyrho (https://github.com/popgenmethods/pyrho)
+        while using population-specific population size history estimates
+        obtained from smc++ (https://github.com/popgenmethods/smcpp).
+        Genetic maps are only available for the 22 autosomes.
+        See https://doi.org/10.1126/sciadv.aaw9206 for more
+        details.""".format(
+            pop
+        ),
+        url=(
+            "https://stdpopsim.s3-us-west-2.amazonaws.com/genetic_maps/"
+            "HomSap/Pyrho{}_GRCh38.tar.gz".format(pop)
+        ),
+        sha256=sha256,
+        file_pattern=("Pyrho{}_GRCh38_chr{{id}}.txt".format(pop)),
+        citations=[
+            stdpopsim.Citation(
+                year=2019,
+                author="Spence and Song",
+                doi="https://doi.org/10.1126/sciadv.aaw9206",
+                reasons={stdpopsim.CiteReason.GEN_MAP},
+            )
+        ],
+    )
+    _species.add_genetic_map(_gm)
diff --git a/stdpopsim/catalog/HomSap/species.py b/stdpopsim/catalog/HomSap/species.py
new file mode 100644
index 000000000..8289f2e19
--- /dev/null
+++ b/stdpopsim/catalog/HomSap/species.py
@@ -0,0 +1,107 @@
+import stdpopsim
+
+from . import genome_data
+
+###########################################################
+#
+# Genome definition
+#
+###########################################################
+
+
+# mean_recombination_rate was computed across all windows of the HapMapII_GRCh37 map
+_recombination_rate_data = {
+    "1": 1.1485597641285933e-08,
+    "2": 1.1054289277533446e-08,
+    "3": 1.1279585624662551e-08,
+    "4": 1.1231162636001008e-08,
+    "5": 1.1280936570022824e-08,
+    "6": 1.1222852661225285e-08,
+    "7": 1.1764614397655721e-08,
+    "8": 1.1478465778920576e-08,
+    "9": 1.1780701596308656e-08,
+    "10": 1.3365134257075317e-08,
+    "11": 1.1719334320833283e-08,
+    "12": 1.305017186986983e-08,
+    "13": 1.0914860554958317e-08,
+    "14": 1.119730771394731e-08,
+    "15": 1.3835785893339787e-08,
+    "16": 1.4834607113882717e-08,
+    "17": 1.582489036239487e-08,
+    "18": 1.5075956950023575e-08,
+    "19": 1.8220141872466202e-08,
+    "20": 1.7178269031631664e-08,
+    "21": 1.3045214034879191e-08,
+    "22": 1.4445022767788226e-08,
+    "X": 1.164662223273842e-08,
+    "Y": 0.0,
+    "MT": 0.0,
+}
+
+_genome2001 = stdpopsim.Citation(
+    doi="http://dx.doi.org/10.1038/35057062",
+    year="2001",
+    author="International Human Genome Sequencing Consortium",
+    reasons={stdpopsim.CiteReason.ASSEMBLY},
+)
+
+_hapmap2007 = stdpopsim.Citation(
+    doi="https://doi.org/10.1038/nature06258",
+    year=2007,
+    author="The International HapMap Consortium",
+)
+
+_takahata1993 = stdpopsim.Citation(
+    doi="https://doi.org/10.1093/oxfordjournals.molbev.a039995",
+    year="1993",
+    author="Takahata",
+    reasons={stdpopsim.CiteReason.POP_SIZE},
+)
+
+_tian2019 = stdpopsim.Citation(
+    doi="https://doi.org/10.1016/j.ajhg.2019.09.012",
+    year="2019",
+    author="Tian, Browning, and Browning",
+    reasons={stdpopsim.CiteReason.MUT_RATE},
+)
+
+_tremblay2000 = stdpopsim.Citation(
+    doi="https://doi.org/10.1086/302770",
+    year="2000",
+    author="Tremblay and Vézina",
+    reasons={stdpopsim.CiteReason.GEN_TIME},
+)
+
+_chromosomes = []
+for name, data in genome_data.data["chromosomes"].items():
+    _chromosomes.append(
+        stdpopsim.Chromosome(
+            id=name,
+            length=data["length"],
+            synonyms=data["synonyms"],
+            mutation_rate=1.29e-8,
+            recombination_rate=_recombination_rate_data[name],
+        )
+    )
+
+_genome = stdpopsim.Genome(
+    chromosomes=_chromosomes,
+    mutation_rate_citations=[_tian2019.because(stdpopsim.CiteReason.MUT_RATE)],
+    recombination_rate_citations=[_hapmap2007.because(stdpopsim.CiteReason.REC_RATE)],
+    assembly_name=genome_data.data["assembly_name"],
+    assembly_accession=genome_data.data["assembly_accession"],
+    assembly_citations=[_genome2001],
+)
+
+_species = stdpopsim.Species(
+    id="HomSap",
+    name="Homo sapiens",
+    common_name="Human",
+    genome=_genome,
+    generation_time=30,
+    generation_time_citations=[_tremblay2000.because(stdpopsim.CiteReason.GEN_TIME)],
+    population_size=10 ** 4,
+    population_size_citations=[_takahata1993.because(stdpopsim.CiteReason.POP_SIZE)],
+)
+
+stdpopsim.register_species(_species)
diff --git a/stdpopsim/catalog/PonAbe/__init__.py b/stdpopsim/catalog/PonAbe/__init__.py
index f4b19c446..fd0977028 100644
--- a/stdpopsim/catalog/PonAbe/__init__.py
+++ b/stdpopsim/catalog/PonAbe/__init__.py
@@ -1,238 +1,6 @@
 """
-Orangutans.
+Catalog data for Pongo abelii
 """
-import math
-
-import msprime
-
-import stdpopsim
-from . import genome_data
-
-###########################################################
-#
-# Genome definition
-#
-###########################################################
-
-# Average recombination rates per chromosome taken from the Pongo abelii
-# recombination map inferred in Nater et al (doi:10.1016/j.cub.2017.09.047)
-# For chromosome X, we use the rate assumed in Locke et al.
-_recombination_rate_data = {
-    "1": 5.19e-9,
-    "2a": 5.43e-9,
-    "2b": 5.38e-9,
-    "3": 5.36e-9,
-    "4": 5.43e-9,
-    "5": 5.20e-9,
-    "6": 5.22e-9,
-    "7": 5.73e-9,
-    "8": 5.67e-9,
-    "9": 5.34e-9,
-    "10": 5.91e-9,
-    "11": 5.29e-9,
-    "12": 5.44e-9,
-    "13": 4.91e-9,
-    "14": 4.70e-9,
-    "15": 4.82e-9,
-    "16": 6.12e-9,
-    "17": 7.26e-9,
-    "18": 4.57e-9,
-    "19": 7.56e-9,
-    "20": 5.83e-9,
-    "21": 4.98e-9,
-    "22": 6.03e-9,
-    "X": 9.50e-9,
-    "MT": 0,
-}
-
-_locke2011 = stdpopsim.Citation(
-    author="Locke et al.", year=2011, doi="http://doi.org/10.1038/nature09687"
-)
-
-_nater2017 = stdpopsim.Citation(
-    author="Nater et al.", year=2017, doi="https://doi.org/10.1016/j.cub.2017.09.047"
-)
-
-_chromosomes = []
-for name, data in genome_data.data["chromosomes"].items():
-    _chromosomes.append(
-        stdpopsim.Chromosome(
-            id=name,
-            length=data["length"],
-            synonyms=data["synonyms"],
-            # Nater et al. 2017 used mu=1.5e-8 per generation, based on the
-            # assumption that it's similar to humans and chimps.
-            mutation_rate=1.5e-8,
-            recombination_rate=_recombination_rate_data[name],
-        )
-    )
-
-_genome = stdpopsim.Genome(
-    chromosomes=_chromosomes,
-    assembly_name=genome_data.data["assembly_name"],
-    assembly_accession=genome_data.data["assembly_accession"],
-    mutation_rate_citations=[_nater2017.because(stdpopsim.CiteReason.MUT_RATE)],
-)
-
-_species = stdpopsim.Species(
-    id="PonAbe",
-    name="Pongo abelii",
-    common_name="Sumatran orangutan",
-    genome=_genome,
-    # generation time used by Locke et al. without further citation
-    generation_time=20,
-    generation_time_citations=[_locke2011.because(stdpopsim.CiteReason.GEN_TIME)],
-    # Locke et al. inferred ancestral Ne
-    population_size=1.79e4,
-    population_size_citations=[_locke2011.because(stdpopsim.CiteReason.POP_SIZE)],
-)
-
-stdpopsim.register_species(_species)
-
-
-###########################################################
-#
-# Genetic maps
-#
-###########################################################
-
-# There are two genetic maps available for Orangutan species: one for Pongo
-# abelii (Sumatran orangutan) and one for Pongo pygmaeus (Bornean orangutan).
-# Both recombination maps were inferred using LDhat in Nater et al. (2017),
-# doi: 10.1016/j.cub.2017.09.047. Both recombination maps are mapped to PonAbe2.
-# Recombination maps from Nater et al. were converted from rho/kbp to cM using
-# Watterson's estimator of theta to estimate Ne = 41,000 (Sumatra) and
-# Ne = 27,000 (Borneo). See supporting information in Nater et al. for details.
-
-
-_gm_pa = stdpopsim.GeneticMap(
-    species=_species,
-    id="NaterPA_PonAbe2",
-    description="From Nater et al. (2017) for Pongo abelii",
-    long_description="""
-        This genetic map is from the Nater et al. (2017) study, inferred using
-        LDhat from n=15 whole-genome sequenced Sumatran orangutan individuals.
-        See https://doi.org/10.1016/j.cub.2017.09.047 for more details.
-        """,
-    url=(
-        "https://stdpopsim.s3-us-west-2.amazonaws.com/genetic_maps/PonAbe/"
-        "NaterPA_PonAbe2.tar.gz"
-    ),
-    sha256="33b0162a6d945dd341bd1086d213ad1bc16949c9210d3b49d92692fd7f831ace",
-    file_pattern="Nater_et_al_PA_chr{id}_PonAbe2.txt",
-    citations=[_nater2017.because(stdpopsim.CiteReason.GEN_MAP)],
-)
-_species.add_genetic_map(_gm_pa)
-
-_gm_pp = stdpopsim.GeneticMap(
-    species=_species,
-    id="NaterPP_PonAbe2",
-    description="From Nater et al. (2017) for Pongo pygmaeus",
-    long_description="""
-        This genetic map is from the Nater et al. (2017) study, inferred using
-        LDhat from n=20 whole-genome sequenced Bornean orangutan individuals.
-        See https://doi.org/10.1016/j.cub.2017.09.047 for more details.
-        """,
-    url=(
-        "https://stdpopsim.s3-us-west-2.amazonaws.com/genetic_maps/PonPyg/"
-        "NaterPP_PonAbe2.tar.gz"
-    ),
-    sha256="f4858b7efe15abe28b9367e7e9dc16614dc614df9326f9eedbcd943ca075bed7",
-    file_pattern="Nater_et_al_PP_chr{id}_PonAbe2.txt",
-    citations=[_nater2017.because(stdpopsim.CiteReason.GEN_MAP)],
-)
-_species.add_genetic_map(_gm_pp)
-
-
-###########################################################
-#
-# Demographic models
-#
-###########################################################
-
-
-def _orangutan():
-    id = "TwoSpecies_2L11"
-    description = "Two population orangutan model"
-    long_description = """
-        The two orang-utan species, Sumatran (Pongo abelii) and Bornean (Pongo
-        pygmaeus) inferred from the joint-site frequency spectrum with ten
-        individuals from each population. This model is an isolation-with-
-        migration model, with exponential growth or decay in each population
-        after the split. The Sumatran population grows in size, while the
-        Bornean population slightly declines.
-    """
-
-    citations = [_locke2011.because(stdpopsim.CiteReason.DEM_MODEL)]
-
-    populations = [
-        stdpopsim.Population("Bornean", "Pongo pygmaeus (Bornean) population"),
-        stdpopsim.Population("Sumatran", "Pongo abelii (Sumatran) population"),
-    ]
-
-    # Parameters from paper:
-    # ancestral size, before split
-    Na = 17934
-
-    # time of split
-    T_split_years = 403149
-    # get split time in units of generations
-    generation_time = _species.generation_time
-    T_split = T_split_years / generation_time
-
-    # proportion of ancestral pop to branch B
-    s = 0.592
-
-    # sizes at split
-    Na_B = 17934 * s
-    Na_S = 17934 * (1 - s)
-
-    # present sizes
-    N_B = 8805
-    N_S = 37661
-
-    # get growth rates
-    r_B = -1 * math.log(Na_B / N_B) / T_split
-    r_S = -1 * math.log(Na_S / N_S) / T_split
-
-    # migration rates
-    m_S_B = 0.395 / 2 / Na
-    m_B_S = 0.239 / 2 / Na
-
-    # Population IDs correspond to their indexes in the population
-    # configuration array. Therefore, we have 0=B and 1=S
-    # initially.
-
-    return stdpopsim.DemographicModel(
-        id=id,
-        description=description,
-        long_description=long_description,
-        citations=citations,
-        populations=populations,
-        generation_time=generation_time,
-        population_configurations=[
-            msprime.PopulationConfiguration(
-                initial_size=N_B, growth_rate=r_B, metadata=populations[0].asdict()
-            ),  # NOQA
-            msprime.PopulationConfiguration(
-                initial_size=N_S, growth_rate=r_S, metadata=populations[1].asdict()
-            ),  # NOQA
-        ],
-        migration_matrix=[
-            [0, m_B_S],  # NOQA
-            [m_S_B, 0],  # NOQA
-        ],
-        demographic_events=[
-            # Merge populations and turn off migration, change to size Na
-            msprime.MassMigration(
-                time=T_split, source=1, destination=0, proportion=1.0
-            ),
-            msprime.MigrationRateChange(time=T_split, rate=0),
-            msprime.PopulationParametersChange(
-                time=T_split, initial_size=Na, growth_rate=0, population_id=0
-            ),
-        ],
-    )
-
-
-_species.add_demographic_model(_orangutan())
+from . import species  # noqa 401
+from . import genetic_maps  # noqa 401
+from . import demographic_models  # noqa 401
diff --git a/stdpopsim/catalog/PonAbe/demographic_models.py b/stdpopsim/catalog/PonAbe/demographic_models.py
new file mode 100644
index 000000000..81f7cd5d3
--- /dev/null
+++ b/stdpopsim/catalog/PonAbe/demographic_models.py
@@ -0,0 +1,100 @@
+import math
+
+import msprime
+
+import stdpopsim
+
+
+_species = stdpopsim.get_species("PonAbe")
+
+
+_locke2011 = stdpopsim.Citation(
+    author="Locke et al.", year=2011, doi="http://doi.org/10.1038/nature09687"
+)
+
+
+def _orangutan():
+    id = "TwoSpecies_2L11"
+    description = "Two population orangutan model"
+    long_description = """
+        The two orang-utan species, Sumatran (Pongo abelii) and Bornean (Pongo
+        pygmaeus) inferred from the joint-site frequency spectrum with ten
+        individuals from each population. This model is an isolation-with-
+        migration model, with exponential growth or decay in each population
+        after the split. The Sumatran population grows in size, while the
+        Bornean population slightly declines.
+    """
+
+    citations = [_locke2011.because(stdpopsim.CiteReason.DEM_MODEL)]
+
+    populations = [
+        stdpopsim.Population("Bornean", "Pongo pygmaeus (Bornean) population"),
+        stdpopsim.Population("Sumatran", "Pongo abelii (Sumatran) population"),
+    ]
+
+    # Parameters from paper:
+    # ancestral size, before split
+    Na = 17934
+
+    # time of split
+    T_split_years = 403149
+    # get split time in units of generations
+    generation_time = _species.generation_time
+    T_split = T_split_years / generation_time
+
+    # proportion of ancestral pop to branch B
+    s = 0.592
+
+    # sizes at split
+    Na_B = 17934 * s
+    Na_S = 17934 * (1 - s)
+
+    # present sizes
+    N_B = 8805
+    N_S = 37661
+
+    # get growth rates
+    r_B = -1 * math.log(Na_B / N_B) / T_split
+    r_S = -1 * math.log(Na_S / N_S) / T_split
+
+    # migration rates
+    m_S_B = 0.395 / 2 / Na
+    m_B_S = 0.239 / 2 / Na
+
+    # Population IDs correspond to their indexes in the population
+    # configuration array. Therefore, we have 0=B and 1=S
+    # initially.
+
+    return stdpopsim.DemographicModel(
+        id=id,
+        description=description,
+        long_description=long_description,
+        citations=citations,
+        populations=populations,
+        generation_time=generation_time,
+        population_configurations=[
+            msprime.PopulationConfiguration(
+                initial_size=N_B, growth_rate=r_B, metadata=populations[0].asdict()
+            ),  # NOQA
+            msprime.PopulationConfiguration(
+                initial_size=N_S, growth_rate=r_S, metadata=populations[1].asdict()
+            ),  # NOQA
+        ],
+        migration_matrix=[
+            [0, m_B_S],  # NOQA
+            [m_S_B, 0],  # NOQA
+        ],
+        demographic_events=[
+            # Merge populations and turn off migration, change to size Na
+            msprime.MassMigration(
+                time=T_split, source=1, destination=0, proportion=1.0
+            ),
+            msprime.MigrationRateChange(time=T_split, rate=0),
+            msprime.PopulationParametersChange(
+                time=T_split, initial_size=Na, growth_rate=0, population_id=0
+            ),
+        ],
+    )
+
+
+_species.add_demographic_model(_orangutan())
diff --git a/stdpopsim/catalog/PonAbe/genetic_maps.py b/stdpopsim/catalog/PonAbe/genetic_maps.py
new file mode 100644
index 000000000..e077e7679
--- /dev/null
+++ b/stdpopsim/catalog/PonAbe/genetic_maps.py
@@ -0,0 +1,54 @@
+import stdpopsim
+
+_species = stdpopsim.get_species("PonAbe")
+
+_nater2017 = stdpopsim.Citation(
+    author="Nater et al.", year=2017, doi="https://doi.org/10.1016/j.cub.2017.09.047"
+)
+
+# There are two genetic maps available for Orangutan species: one for Pongo
+# abelii (Sumatran orangutan) and one for Pongo pygmaeus (Bornean orangutan).
+# Both recombination maps were inferred using LDhat in Nater et al. (2017),
+# doi: 10.1016/j.cub.2017.09.047. Both recombination maps are mapped to PonAbe2.
+# Recombination maps from Nater et al. were converted from rho/kbp to cM using
+# Watterson's estimator of theta to estimate Ne = 41,000 (Sumatra) and
+# Ne = 27,000 (Borneo). See supporting information in Nater et al. for details.
+
+
+_gm_pa = stdpopsim.GeneticMap(
+    species=_species,
+    id="NaterPA_PonAbe2",
+    description="From Nater et al. (2017) for Pongo abelii",
+    long_description="""
+        This genetic map is from the Nater et al. (2017) study, inferred using
+        LDhat from n=15 whole-genome sequenced Sumatran orangutan individuals.
+        See https://doi.org/10.1016/j.cub.2017.09.047 for more details.
+        """,
+    url=(
+        "https://stdpopsim.s3-us-west-2.amazonaws.com/genetic_maps/PonAbe/"
+        "NaterPA_PonAbe2.tar.gz"
+    ),
+    sha256="33b0162a6d945dd341bd1086d213ad1bc16949c9210d3b49d92692fd7f831ace",
+    file_pattern="Nater_et_al_PA_chr{id}_PonAbe2.txt",
+    citations=[_nater2017.because(stdpopsim.CiteReason.GEN_MAP)],
+)
+_species.add_genetic_map(_gm_pa)
+
+_gm_pp = stdpopsim.GeneticMap(
+    species=_species,
+    id="NaterPP_PonAbe2",
+    description="From Nater et al. (2017) for Pongo pygmaeus",
+    long_description="""
+        This genetic map is from the Nater et al. (2017) study, inferred using
+        LDhat from n=20 whole-genome sequenced Bornean orangutan individuals.
+        See https://doi.org/10.1016/j.cub.2017.09.047 for more details.
+        """,
+    url=(
+        "https://stdpopsim.s3-us-west-2.amazonaws.com/genetic_maps/PonPyg/"
+        "NaterPP_PonAbe2.tar.gz"
+    ),
+    sha256="f4858b7efe15abe28b9367e7e9dc16614dc614df9326f9eedbcd943ca075bed7",
+    file_pattern="Nater_et_al_PP_chr{id}_PonAbe2.txt",
+    citations=[_nater2017.because(stdpopsim.CiteReason.GEN_MAP)],
+)
+_species.add_genetic_map(_gm_pp)
diff --git a/stdpopsim/catalog/PonAbe/species.py b/stdpopsim/catalog/PonAbe/species.py
new file mode 100644
index 000000000..3c80bc528
--- /dev/null
+++ b/stdpopsim/catalog/PonAbe/species.py
@@ -0,0 +1,78 @@
+import stdpopsim
+from . import genome_data
+
+
+# Average recombination rates per chromosome taken from the Pongo abelii
+# recombination map inferred in Nater et al (doi:10.1016/j.cub.2017.09.047)
+# For chromosome X, we use the rate assumed in Locke et al.
+_recombination_rate_data = {
+    "1": 5.19e-9,
+    "2a": 5.43e-9,
+    "2b": 5.38e-9,
+    "3": 5.36e-9,
+    "4": 5.43e-9,
+    "5": 5.20e-9,
+    "6": 5.22e-9,
+    "7": 5.73e-9,
+    "8": 5.67e-9,
+    "9": 5.34e-9,
+    "10": 5.91e-9,
+    "11": 5.29e-9,
+    "12": 5.44e-9,
+    "13": 4.91e-9,
+    "14": 4.70e-9,
+    "15": 4.82e-9,
+    "16": 6.12e-9,
+    "17": 7.26e-9,
+    "18": 4.57e-9,
+    "19": 7.56e-9,
+    "20": 5.83e-9,
+    "21": 4.98e-9,
+    "22": 6.03e-9,
+    "X": 9.50e-9,
+    "MT": 0,
+}
+
+_locke2011 = stdpopsim.Citation(
+    author="Locke et al.", year=2011, doi="http://doi.org/10.1038/nature09687"
+)
+
+_nater2017 = stdpopsim.Citation(
+    author="Nater et al.", year=2017, doi="https://doi.org/10.1016/j.cub.2017.09.047"
+)
+
+_chromosomes = []
+for name, data in genome_data.data["chromosomes"].items():
+    _chromosomes.append(
+        stdpopsim.Chromosome(
+            id=name,
+            length=data["length"],
+            synonyms=data["synonyms"],
+            # Nater et al. 2017 used mu=1.5e-8 per generation, based on the
+            # assumption that it's similar to humans and chimps.
+            mutation_rate=1.5e-8,
+            recombination_rate=_recombination_rate_data[name],
+        )
+    )
+
+_genome = stdpopsim.Genome(
+    chromosomes=_chromosomes,
+    assembly_name=genome_data.data["assembly_name"],
+    assembly_accession=genome_data.data["assembly_accession"],
+    mutation_rate_citations=[_nater2017.because(stdpopsim.CiteReason.MUT_RATE)],
+)
+
+_species = stdpopsim.Species(
+    id="PonAbe",
+    name="Pongo abelii",
+    common_name="Sumatran orangutan",
+    genome=_genome,
+    # generation time used by Locke et al. without further citation
+    generation_time=20,
+    generation_time_citations=[_locke2011.because(stdpopsim.CiteReason.GEN_TIME)],
+    # Locke et al. inferred ancestral Ne
+    population_size=1.79e4,
+    population_size_citations=[_locke2011.because(stdpopsim.CiteReason.POP_SIZE)],
+)
+
+stdpopsim.register_species(_species)


Test patch:
diff --git a/tests/test_AraTha.py b/tests/test_AraTha.py
index 00780b895..40400f66f 100644
--- a/tests/test_AraTha.py
+++ b/tests/test_AraTha.py
@@ -1,18 +1,16 @@
-import unittest
-
 import stdpopsim
 from tests import test_species
 
 
-class TestSpecies(unittest.TestCase, test_species.SpeciesTestMixin):
+class TestSpecies(test_species.SpeciesTestBase):
     species = stdpopsim.get_species("AraTha")
 
     def test_basic_attributes(self):
-        self.assertEqual(self.species.population_size, 10 ** 4)
-        self.assertEqual(self.species.generation_time, 1)
+        assert self.species.population_size == 10 ** 4
+        assert self.species.generation_time == 1
 
 
-class TestGenome(unittest.TestCase, test_species.GenomeTestMixin):
+class TestGenome(test_species.GenomeTestBase):
     """
     Tests for the arabidopsis_thaliana genome.
     """
@@ -21,4 +19,4 @@ class TestGenome(unittest.TestCase, test_species.GenomeTestMixin):
 
     def test_basic_attributes(self):
         # 5 autosomes + Mt + Pt
-        self.assertEqual(len(self.genome.chromosomes), 7)
+        assert len(self.genome.chromosomes) == 7
diff --git a/tests/test_CanFam.py b/tests/test_CanFam.py
index 38d8de1d3..45a09336c 100644
--- a/tests/test_CanFam.py
+++ b/tests/test_CanFam.py
@@ -1,31 +1,18 @@
-import os
-import unittest
-from unittest import mock
-
 import stdpopsim
-import stdpopsim.cli
 from tests import test_species
 
 
-class TestSpecies(unittest.TestCase, test_species.SpeciesTestMixin):
+class TestSpecies(test_species.SpeciesTestBase):
     species = stdpopsim.get_species("CanFam")
 
     def test_basic_attributes(self):
-        self.assertEqual(self.species.population_size, 13000)
-        self.assertEqual(self.species.generation_time, 3)
+        assert self.species.population_size == 13000
+        assert self.species.generation_time == 3
 
 
-class TestGenome(unittest.TestCase, test_species.GenomeTestMixin):
+class TestGenome(test_species.GenomeTestBase):
     genome = stdpopsim.get_species("CanFam").genome
 
     def test_basic_attributes(self):
         nchrom = 40  # 38 + X + MT
-        self.assertEqual(len(self.genome.chromosomes), nchrom)
-
-
-class TestThatDogsCanBeSimulated(unittest.TestCase):
-    def test_basic_cli_usage(self):
-        cmd = "-q CanFam -c 38 -l 0.001 --seed 1234 10"
-        with mock.patch("sys.stdout", autospec=True) as stdout:
-            stdout.buffer = open(os.devnull, "wb")
-            stdpopsim.cli.stdpopsim_main(cmd.split())
+        assert len(self.genome.chromosomes) == nchrom
diff --git a/tests/test_DroMel.py b/tests/test_DroMel.py
index 4878d34c9..66718e8c3 100644
--- a/tests/test_DroMel.py
+++ b/tests/test_DroMel.py
@@ -1,21 +1,21 @@
 """
 Tests for the drosophila_melanogaster data definitions.
 """
-import unittest
+import pytest
 
 import stdpopsim
 from tests import test_species
 
 
-class TestSpecies(unittest.TestCase, test_species.SpeciesTestMixin):
+class TestSpecies(test_species.SpeciesTestBase):
     species = stdpopsim.get_species("DroMel")
 
     def test_basic_attributes(self):
-        self.assertEqual(self.species.population_size, 1720600)
-        self.assertEqual(self.species.generation_time, 0.1)
+        self.species.population_size == 1720600
+        self.species.generation_time == 0.1
 
 
-class TestGenome(unittest.TestCase, test_species.GenomeTestMixin):
+class TestGenome(test_species.GenomeTestBase):
     """
     Tests for the drosophila_melanogaster genome.
     """
@@ -23,15 +23,14 @@ class TestGenome(unittest.TestCase, test_species.GenomeTestMixin):
     genome = stdpopsim.get_species("DroMel").genome
 
     def test_basic_attributes(self):
-        self.assertEqual(len(self.genome.chromosomes), 8)
-
-    def test_non_recombining_chrs(self):
-        chrom = self.genome.get_chromosome("Y")
-        self.assertEqual(chrom.recombination_rate, 0)
-        chrom = self.genome.get_chromosome("mitochondrion_genome")
-        self.assertEqual(chrom.recombination_rate, 0)
-
-    def test_recombining_chrs(self):
-        for name in ["2L", "2R", "3L", "3R", "4", "X"]:
-            chrom = self.genome.get_chromosome(name)
-            self.assertGreater(chrom.recombination_rate, 0)
+        assert len(self.genome.chromosomes) == 8
+
+    @pytest.mark.parametrize("chr_id", ["Y", "mitochondrion_genome"])
+    def test_non_recombining_chrs(self, chr_id):
+        chrom = self.genome.get_chromosome(chr_id)
+        assert chrom.recombination_rate == 0
+
+    @pytest.mark.parametrize("chr_id", ["2L", "2R", "3L", "3R", "4", "X"])
+    def test_recombining_chrs(self, chr_id):
+        chrom = self.genome.get_chromosome(chr_id)
+        assert chrom.recombination_rate > 0
diff --git a/tests/test_EscCol.py b/tests/test_EscCol.py
index 764867bf2..4657a8fcb 100644
--- a/tests/test_EscCol.py
+++ b/tests/test_EscCol.py
@@ -1,25 +1,23 @@
 """
 Tests for the e. coli data definitions.
 """
-import unittest
-
 import stdpopsim
 from tests import test_species
 
 
-class TestSpecies(unittest.TestCase, test_species.SpeciesTestMixin):
+class TestSpecies(test_species.SpeciesTestBase):
     species = stdpopsim.get_species("EscCol")
 
     def test_basic_attributes(self):
         # From paper https://doi.org/10.1093/molbev/msw048
         # Ne taken from Table 2
-        self.assertEqual(self.species.population_size, 1.8e8)
+        assert self.species.population_size == 1.8e8
         # 20 minutes per generation
         generation_time = 1.0 / (525600 / 20)
-        self.assertAlmostEqual(self.species.generation_time, generation_time)
+        assert round(abs(self.species.generation_time - generation_time), 7) == 0
 
 
-class TestGenome(unittest.TestCase, test_species.GenomeTestMixin):
+class TestGenome(test_species.GenomeTestBase):
     """
     Tests for the e_coli genome.
     """
@@ -27,4 +25,4 @@ class TestGenome(unittest.TestCase, test_species.GenomeTestMixin):
     genome = stdpopsim.get_species("EscCol").genome
 
     def test_basic_attributes(self):
-        self.assertEqual(len(self.genome.chromosomes), 1)
+        assert len(self.genome.chromosomes) == 1
diff --git a/tests/test_HomSap.py b/tests/test_HomSap.py
index 131826f97..2e769bf8b 100644
--- a/tests/test_HomSap.py
+++ b/tests/test_HomSap.py
@@ -1,44 +1,48 @@
 """
 Tests for the human data definitions.
 """
-import unittest
-
 import stdpopsim
 
+import pytest
+
 from tests import test_species
 
 
-class TestSpecies(unittest.TestCase, test_species.SpeciesTestMixin):
+class TestSpecies(test_species.SpeciesTestBase):
     species = stdpopsim.get_species("HomSap")
 
     def test_basic_attributes(self):
-        self.assertEqual(self.species.population_size, 10 ** 4)
-        self.assertEqual(self.species.generation_time, 30)
+        assert self.species.population_size == 10 ** 4
+        assert self.species.generation_time == 30
 
 
-class TestGenome(unittest.TestCase, test_species.GenomeTestMixin):
-    """
-    Tests for the human genome.
-    """
+class TestGenome(test_species.GenomeTestBase):
 
     genome = stdpopsim.get_species("HomSap").genome
 
     def test_basic_attributes(self):
-        self.assertEqual(len(self.genome.chromosomes), 25)
+        assert len(self.genome.chromosomes) == 25
 
-    def test_recombination_rates(self):
+    @pytest.mark.parametrize("chr_id", [chrom.id for chrom in genome.chromosomes])
+    def test_recombination_rates(self, chr_id):
         # recompute recombination rates from HapMapII_GRCh37 map then
         # compare the results to the current recombination rates for each chromosome
         genetic_map = "HapMapII_GRCh37"
         species = stdpopsim.get_species("HomSap")
-        for chrom in self.genome.chromosomes:
-            if chrom.id == "Y":
-                with self.assertWarns(Warning):
-                    contig = species.get_contig(chrom.id, genetic_map=genetic_map)
-                print("HERE")
-            else:
-                contig = species.get_contig(chrom.id, genetic_map=genetic_map)
-            self.assertAlmostEqual(
-                chrom.recombination_rate,
-                contig.recombination_map.mean_rate,
-            )
+        chrom = species.genome.get_chromosome(chr_id)
+        if chr_id in ["X", "Y", "MT"]:
+            with pytest.warns(stdpopsim.NonAutosomalWarning):
+                contig = species.get_contig(chr_id, genetic_map=genetic_map)
+        elif chr_id in ["3", "5", "7", "11", "16", "17", "18", "20"]:
+            contig = species.get_contig(chr_id, genetic_map=genetic_map)
+        else:
+            # The rest of the chromosomes are currently emitting a warning about
+            # the mismatch in chromosome lengths because of the fact that we're
+            # on 37 for the map. This should be resolved when we start using the
+            # lifted over map.
+            with pytest.warns(UserWarning, match="longer than chromosome length"):
+                contig = species.get_contig(chr_id, genetic_map=genetic_map)
+        assert pytest.approx(
+            chrom.recombination_rate,
+            contig.recombination_map.mean_rate,
+        )
diff --git a/tests/test_PonAbe.py b/tests/test_PonAbe.py
index 86aded10a..c04990bdb 100644
--- a/tests/test_PonAbe.py
+++ b/tests/test_PonAbe.py
@@ -1,18 +1,16 @@
-import unittest
-
 import stdpopsim
 from tests import test_species
 
 
-class TestSpecies(unittest.TestCase, test_species.SpeciesTestMixin):
+class TestSpecies(test_species.SpeciesTestBase):
     species = stdpopsim.get_species("PonAbe")
 
     def test_basic_attributes(self):
-        self.assertEqual(self.species.population_size, 1.79 * 10 ** 4)
-        self.assertEqual(self.species.generation_time, 20)
+        assert self.species.population_size == 1.79 * 10 ** 4
+        assert self.species.generation_time == 20
 
 
-class TestGenome(unittest.TestCase, test_species.GenomeTestMixin):
+class TestGenome(test_species.GenomeTestBase):
     """
     Tests for the Pongo abelii genome.
     """
@@ -21,4 +19,4 @@ class TestGenome(unittest.TestCase, test_species.GenomeTestMixin):
 
     def test_basic_attributes(self):
         # 23 + X + MT
-        self.assertEqual(len(self.genome.chromosomes), 25)
+        assert len(self.genome.chromosomes) == 25
diff --git a/tests/test_cli.py b/tests/test_cli.py
index 678849163..a4d4eee6f 100644
--- a/tests/test_cli.py
+++ b/tests/test_cli.py
@@ -16,12 +16,13 @@
 import tskit
 import msprime
 import kastore
+import pytest
 
 import stdpopsim
 import stdpopsim.cli as cli
 
 
-class TestException(Exception):
+class ExceptionForTesting(Exception):
     """
     Custom exception we can throw for testing.
     """
@@ -513,9 +514,9 @@ def run_stdpopsim(self, command, mock_setup_logging):
 
     def test_exit(self):
         with mock.patch(
-            "sys.exit", side_effect=TestException, autospec=True
+            "sys.exit", side_effect=ExceptionForTesting, autospec=True
         ) as mocked_exit:
-            with self.assertRaises(TestException):
+            with self.assertRaises(ExceptionForTesting):
                 cli.exit("XXX")
             mocked_exit.assert_called_once()
             args = mocked_exit.call_args[0]
@@ -527,9 +528,9 @@ def test_exit(self):
     @mock.patch("stdpopsim.cli.setup_logging", autospec=True)
     def verify_bad_samples(self, cmd, mock_setup_logging):
         with mock.patch(
-            "stdpopsim.cli.exit", side_effect=TestException, autospec=True
+            "stdpopsim.cli.exit", side_effect=ExceptionForTesting, autospec=True
         ) as mocked_exit:
-            with self.assertRaises(TestException):
+            with self.assertRaises(ExceptionForTesting):
                 cli.stdpopsim_main(cmd.split())
             mocked_exit.assert_called_once()
 
@@ -549,9 +550,11 @@ def test_browning_america(self):
 class TestHelp(unittest.TestCase):
     def run_stdpopsim(self, command):
         with mock.patch(
-            "argparse.ArgumentParser.exit", side_effect=TestException, autospec=True
+            "argparse.ArgumentParser.exit",
+            side_effect=ExceptionForTesting,
+            autospec=True,
         ) as mocked_exit:
-            with self.assertRaises(TestException):
+            with self.assertRaises(ExceptionForTesting):
                 capture_output(cli.stdpopsim_main, command.split())
             mocked_exit.assert_called_once()
 
@@ -931,3 +934,14 @@ def test_noQC_citations_not_written(self):
 
     def test_noQC_citations_not_written_verbose(self):
         self.verify_noQC_citations_not_written("-vv EscCol -d FakeModel -D 10 -L 10")
+
+
+@pytest.mark.parametrize(
+    "species_id", [species.id for species in stdpopsim.all_species()]
+)
+def test_species_simulation(species_id):
+    cmd = f"-q {species_id} -L 1 --seed 1234 10"
+    # Just check to see if the simulation runs
+    with mock.patch("sys.stdout", autospec=True) as stdout:
+        stdout.buffer = open(os.devnull, "wb")
+        stdpopsim.cli.stdpopsim_main(cmd.split())
diff --git a/tests/test_species.py b/tests/test_species.py
index 604bbdbec..7b3e246aa 100644
--- a/tests/test_species.py
+++ b/tests/test_species.py
@@ -1,5 +1,5 @@
 """
-Tests for the genetic species interface.
+Tests for the generic species interface.
 """
 import unittest
 import math
@@ -88,42 +88,42 @@ def test_add_duplicate_annotation(self):
             species.add_annotations(an)
 
 
-class SpeciesTestMixin(object):
+class SpeciesTestBase:
     """
-    Mixin class for testing individual species properties.
+    Base class for testing individual species properties.
     """
 
     species = None  # To be defined in subclasses.
 
     def test_str(self):
         s = str(self.species)
-        self.assertGreater(len(s), 0)
-        self.assertIsInstance(s, str)
+        assert len(s) > 0
+        assert isinstance(s, str)
 
     def test_id(self):
-        self.assertIsInstance(self.species.id, str)
-        self.assertTrue(utils.is_valid_species_id(self.species.id))
+        assert isinstance(self.species.id, str)
+        assert utils.is_valid_species_id(self.species.id)
 
     def test_name(self):
-        self.assertIsInstance(self.species.name, str)
-        self.assertTrue(utils.is_valid_species_name(self.species.name))
+        assert isinstance(self.species.name, str)
+        assert utils.is_valid_species_name(self.species.name)
 
     def test_common_name(self):
-        self.assertIsInstance(self.species.name, str)
-        self.assertTrue(utils.is_valid_species_common_name(self.species.common_name))
+        assert isinstance(self.species.common_name, str)
+        assert utils.is_valid_species_common_name(self.species.name)
 
 
-class GenomeTestMixin(object):
+class GenomeTestBase:
     """
-    Mixin class for testing individual genome properties.
+    Base class for testing individual genome properties.
     """
 
     genome = None  # To be defined in subclasses.
 
     def test_str(self):
         s = str(self.genome)
-        self.assertGreater(len(s), 0)
-        self.assertIsInstance(s, str)
+        assert len(s) > 0
+        assert isinstance(s, str)
 
     def test_mean_recombination_rate(self):
         # test that the mean recombination rate lies between the max and min values
@@ -137,8 +137,8 @@ def test_mean_recombination_rate(self):
             highest_rr = max(highest_rr, rr)
         mean_genome_rr = self.genome.mean_recombination_rate
         if not math.isclose(mean_genome_rr, lowest_rr):
-            self.assertGreaterEqual(mean_genome_rr, lowest_rr)
-            self.assertGreaterEqual(highest_rr, mean_genome_rr)
+            assert mean_genome_rr >= lowest_rr
+            assert highest_rr >= mean_genome_rr
 
     def test_mean_mutation_rate(self):
         # test that the mean mutation rate lies between the max and min values
@@ -150,8 +150,8 @@ def test_mean_mutation_rate(self):
             highest_mr = max(highest_mr, mr)
         mean_genome_mr = self.genome.mean_mutation_rate
         if not math.isclose(mean_genome_mr, lowest_mr):
-            self.assertGreaterEqual(mean_genome_mr, lowest_mr)
-            self.assertGreaterEqual(highest_mr, mean_genome_mr)
+            assert mean_genome_mr >= lowest_mr
+            assert highest_mr >= mean_genome_mr
 
 
 class TestAllGenomes(unittest.TestCase):