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

Issue:
What to do when our sims don't line up with data?
Okay so just ran into an interesting case with @cjbattey.

CJ has been comparing stdpopsim to 1000 Genomes data as a baseline for a project we are working on. Simulating under `OutOfAfrica_3G09` yielded too few mutations in comparison to equivalent data from 1kg. e.g., for chr22 and considering only sites in the 'P' state of the 1kg masking `OutOfAfrica_3G09` gives roughly 2x fewer SNPs than are seen in phase 3 samples.

One reason for this difference is that stdpopsim default mutation rate (1.29e-08) is almost half of the mutation rate estimate used by @RyanGutenkunst et al in the original paper (2.38e-08). 

So this is an issue that we might deal with in a number of ways (see #530).

1. have paper specific mutation rates override defaults (as others have suggested)
2. trigger warnings if defaults are used (our rationale for having defaults is still good-- we want up to date, best info in those slots, but users need to beware of caveats)
3. allow users to simply define there own custom parameters. perhaps we would want this to extend to the CLI. i just submitted #556 as part of this.

thoughts? it's kinda thorny

unfroze attributes in engine, contig, chromosome
For a project with @cjbattey we were messing with mutation rates to try to get stdpopsim models to lineup with empirical data. This is a whole other can of worms, but in the meanwhile I unfroze the objects referred to in #530 


What to do when our sims don't line up with data?
Okay so just ran into an interesting case with @cjbattey.

CJ has been comparing stdpopsim to 1000 Genomes data as a baseline for a project we are working on. Simulating under `OutOfAfrica_3G09` yielded too few mutations in comparison to equivalent data from 1kg. e.g., for chr22 and considering only sites in the 'P' state of the 1kg masking `OutOfAfrica_3G09` gives roughly 2x fewer SNPs than are seen in phase 3 samples.

One reason for this difference is that stdpopsim default mutation rate (1.29e-08) is almost half of the mutation rate estimate used by @RyanGutenkunst et al in the original paper (2.38e-08). 

So this is an issue that we might deal with in a number of ways (see #530).

1. have paper specific mutation rates override defaults (as others have suggested)
2. trigger warnings if defaults are used (our rationale for having defaults is still good-- we want up to date, best info in those slots, but users need to beware of caveats)
3. allow users to simply define there own custom parameters. perhaps we would want this to extend to the CLI. i just submitted #556 as part of this.

thoughts? it's kinda thorny



Discussion / hints:
I think users should probably use the Python API and override the values at that point. I don't think we can accommodate that level of detail in the CLI.
The other option is

4. Have "mutation rate used for inference" be part of the model, and then rescale times accordingly.

This is thorny because in some methods the mutation rate isn't just a simple time scaling, though. 
# [Codecov](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/556?src=pr&el=h1) Report
> Merging [#556](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/556?src=pr&el=desc) into [master](https://codecov.io/gh/popsim-consortium/stdpopsim/commit/f6c583d7fb04a3f05afcd8aab963a28bc827e0b6&el=desc) will **not change** coverage.
> The diff coverage is `100.00%`.

[![Impacted file tree graph](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/556/graphs/tree.svg?width=650&height=150&src=pr&token=AOkcaVfvhX)](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/556?src=pr&el=tree)

```diff
@@           Coverage Diff           @@
##           master     #556   +/-   ##
=======================================
  Coverage   99.69%   99.69%           
=======================================
  Files          22       22           
  Lines        1949     1949           
  Branches      202      202           
=======================================
  Hits         1943     1943           
  Misses          2        2           
  Partials        4        4           
```


| [Impacted Files](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/556?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [stdpopsim/engines.py](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/556/diff?src=pr&el=tree#diff-c3RkcG9wc2ltL2VuZ2luZXMucHk=) | `100.00% <100.00%> (ø)` | |
| [stdpopsim/genomes.py](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/556/diff?src=pr&el=tree#diff-c3RkcG9wc2ltL2dlbm9tZXMucHk=) | `100.00% <100.00%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/556?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/556?src=pr&el=footer). Last update [f6c583d...f6757fa](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/556?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

Thanks @andrewkern. frozen=False is actually the default in the attrs library, so you can just specify the decorator as `@attr.s` (without any brackets) and get the same effect. Do you mind changing the diff to have that instead?
how's this?
Perfect, thanks!
I think users should probably use the Python API and override the values at that point. I don't think we can accommodate that level of detail in the CLI.
The other option is

4. Have "mutation rate used for inference" be part of the model, and then rescale times accordingly.

This is thorny because in some methods the mutation rate isn't just a simple time scaling, though. 

Fix patch:
diff --git a/docs/development.rst b/docs/development.rst
index cc52f0557..14503f734 100644
--- a/docs/development.rst
+++ b/docs/development.rst
@@ -563,6 +563,7 @@ The demographic model function should follow this format:
         ]
 
         generation_time = "FILL ME"
+        mutation_rate = "FILL ME"
 
         # parameter value definitions based on published values
 
@@ -573,6 +574,7 @@ The demographic model function should follow this format:
             populations=populations,
             citations=citations,
             generation_time=generation_time,
+            mutation_rate=mutation_rate,
             population_configurations=["FILL ME"],
             migration_matrix=["FILL ME"],
             demographic_events=["FILL ME"],
@@ -602,7 +604,10 @@ The demographic model should include the following:
   a specified reason for citing this model.
 * ``generation_time``: The generation time for the species in years. If you are
   implementing a generic model, the generation time should default to 1.
-
+* ``mutation_rate``: The mutation rate assumed during the inference of this demographic
+  model, if a mutation rate was used. If no mutation rate is associated with this
+  demographic model, which is generally uncommon but possible, depending on the
+  inference method, the mutation rate should be set to ``None``.
 
 Every demographic model has a few necessary features or attributes. First of all,
 demographic models are defined by the population sizes, migration rates, split and
@@ -664,6 +669,8 @@ Now check that your new demographic model function has been imported:
     # Zigzag_1S14
     # AncientEurasia_9K19
     # PapuansOutOfAfrica_10J19
+    # AshkSub_7G19
+    # OutOfAfrica_4J17
 
 
 The example above lists the imported demographic models for humans.
diff --git a/docs/tutorial.rst b/docs/tutorial.rst
index 39dbfa6fa..a0643c999 100644
--- a/docs/tutorial.rst
+++ b/docs/tutorial.rst
@@ -75,7 +75,9 @@ chromosome, recombination map, and number of samples.
 The most basic simulation we can run is to simulate two (haploid) genomes
 - i.e., two samples -
 using the species' defaults as seen in the species help (``stdpopsim HomSap --help``).
-These defaults include constant size population and a uniform recombination map.
+These defaults include constant size population, a uniform recombination map based
+on the average recombination rate (either genome-wide or within a chromosome, if
+specified), and the mutation rate shown above.
 To save time we will specify that the simulation use
 chromosome 22, using the ``-c`` option. We also specify that the resulting
 tree-sequence formatted output should be written to the file ``foo.ts`` with the
@@ -124,6 +126,15 @@ The order of those numbers is the same as the order
 specified in the model documentation. In this case, ``2 3`` means
 we are simulating two African American samples and three European American samples.
 
+.. note::
+    Many demographic models were inferred or calibrated using a mutation rate that
+    differs from the cataloged species' mutation rate. Simulations using the CLI now
+    automatically use the *model's* specified mutation rate instead of the species
+    rate, so that expected levels of diversity more closely match those observed in
+    the data that were used to infer the demographic model. For generic demographic
+    models or those without associated mutation rates, the species mutation rate is
+    used.
+
 Now we want to add an empirical recombination map to make the simulation more
 realistic. We can look up the available recombination maps using the
 ``--help-genetic-maps`` flag (here, truncated for space):
@@ -282,7 +293,7 @@ using SLiM with a (very large!) scaling factor of 1000, we could run
 
 .. code-block:: console
 
-   $ stdpopsim -e slim --slim-scaling-factor 1000 DroMel \
+   $ stdpopsim -e slim --slim-scaling-factor 1000 DroMel\
    $     -c chr2L -l 0.05 -o foo.ts -d African3Epoch_1S16 100
 
 The scaling factor of 1000 makes this model run very quickly,
@@ -350,31 +361,21 @@ Running a published model
 The first example uses a mostly default genome
 with a published demographic model.
 
-Pick a species and a contig
----------------------------
+Pick a species and demographic model
+------------------------------------
 
-First, we will pick a species (humans), contig (chromosome 22),
-and genetic map (a flat map, the default):
+First, we will pick a species (here, humans) and the published demographic
+model to simulated under. In ``stdpopsim`` there are two types of model: ones
+taken to match the :ref:`demographic history reported in published papers
+<sec_catalog>`, and :ref:`"generic" models <sec_api_generic_models>`. We'll
+first simulate using a published model from the catalog. Let's see what
+demographic models are available for humans:
 
 .. code-block:: python
 
    import stdpopsim
 
    species = stdpopsim.get_species("HomSap")
-   contig = species.get_contig("chr22")  # default is a flat genetic map
-
-Choose a demographic model
---------------------------
-
-Now, we need a demographic model.
-In ``stdpopsim`` there are two types of model:
-ones taken to match the :ref:`demographic history reported in published papers <sec_catalog>`,
-and :ref:`"generic" models <sec_api_generic_models>`.
-First, we'll simulate samples from a published model of human prehistory.
-Let's see what demographic models are available
-for humans:
-
-.. code-block:: python
 
    for x in species.demographic_models:
        print(x.id)
@@ -387,11 +388,13 @@ for humans:
    # Zigzag_1S14
    # AncientEurasia_9K19
    # PapuansOutOfAfrica_10J19
+   # AshkSub_7G19
+   # OutOfAfrica_4J17
 
 These models are described in detail in the :ref:`Catalog <sec_catalog>`.
 We'll look at the first model, "OutOfAfrica_3G09", from
 `Gutenkunst et al (2009) <https://doi.org/10.1371/journal.pgen.1000695>`__.
-First we need to know how many populations it has, and what they are:
+We can check how many populations exist in this model, and what they are:
 
 .. code-block:: python
 
@@ -400,15 +403,51 @@ First we need to know how many populations it has, and what they are:
    # 3
    print(model.num_sampling_populations)
    # 3
-   print([pop.id for pop in model.populations])
+   print([pop.name for pop in model.populations])
    # ['YRI', 'CEU', 'CHB']
 
-This model has 3 populations, named YRI, CEU, CHB,
-and all three can be sampled from.
+This model has 3 populations, named YRI, CEU, CHB, and all three can be sampled from.
 The number of "sampling" populations could be smaller than the number of populations,
 since some models have ancient populations which are currently not allowed to be
 sampled from - but that is not the case in this model.
 
+Set up the contig
+-----------------
+
+We'll next define the contig, which contains information about the genome length we
+want to simulate and recombination and mutation rates. Here, we use the human
+chromosome 22. If no recombination map is specified, we assume a uniform genetic map
+based on the average recombination rate for that chromosome.
+
+.. code-block:: python
+
+   contig = species.get_contig("chr22")
+
+   # default is a flat genetic map
+   print("mean recombination rate:", f"{contig.recombination_map.mean_rate:.3}")
+   # mean recombination rate: 1.44e-08
+
+   # and the default mutation rate is based on the species default
+   print("mean mutation rate:", contig.mutation_rate)
+   # mean mutation rate: 1.29e-08
+
+   # but note that the mutation rate differs from the model's assumed rate
+   print("model mutation rate:", model.mutation_rate)
+   # model mutation rate: 2.35e-08
+
+The Gutenkunst OOA model was inferred using a mutation rate much larger than the
+default mutation rate in the `stdpopsim` catalog. As such, simulating using this
+model and default rate will result in levels of diversity substantially lower than
+expected for the human population data that this model was inferred from. To match
+observed diversity in humans, we should instead use the mutation rate associated
+with the demographic model:
+
+.. code-block:: python
+
+   contig = species.get_contig("chr22", mutation_rate=model.mutation_rate)
+   print(contig.mutation_rate == model.mutation_rate)
+   # True
+
 Choose a sampling scheme and simulate
 -------------------------------------
 
@@ -459,19 +498,16 @@ so we use the json module to easily parse it:
 Running a generic model
 =======================
 
-Next, we will simulate using a "generic" model,
-with piecewise constant population size.
-This time,
-we will simulate a small portion of human chromosome 2,
-again with a flat recombination map,
-using the current best estimate of the human
-effective population size from the :ref:`sec_catalog`.
+Next, we will simulate using a "generic" model, with piecewise constant
+population size. This time, we will simulate a given genome length under
+a flat recombination map, using an estimate of the
+human effective population size from the :ref:`sec_catalog`.
 
 Choose a species
 ----------------
 
-Although the model is generic, we still need a species,
-to get contig information.
+Although the model is generic, we still need a species in order
+to get the contig information.
 Again, we'll use `Homo sapiens`, which has the id "HomSap".
 (But, you could use any species from the :ref:`sec_catalog`!)
 
@@ -481,25 +517,6 @@ Again, we'll use `Homo sapiens`, which has the id "HomSap".
 
     species = stdpopsim.get_species("HomSap")
 
-Choose a contig and recombination map
--------------------------------------
-
-Next, we set the contig information, to 5% of chromosome 2,
-about 12Mb. Again, you could use a fraction of any of the
-chromosomes listed in the :ref:`sec_catalog`, keeping in mind that
-larger contigs will take longer to simulate.
-
-.. code-block:: python
-
-    contig = species.get_contig("chr2", length_multiplier=0.05)
-    print(contig.recombination_map.get_sequence_length())
-    # 12159968.65
-
-The "sequence length" is the length in base pairs.
-It is not an integer because msprime works in continuous genomic coordinates.
-Again, since we didn't specify a genetic map, by default
-a "flat" map of constant recombination rate is used.
-
 Set up the generic model
 ------------------------
 
@@ -515,6 +532,32 @@ single population of constant over all time.
 Each species has a "default" population size, ``species.population_size``,
 which for humans is 10,000.
 
+Choose a contig and recombination map
+-------------------------------------
+
+Next, we set the contig information. Again, we could use any of the chromosomes
+listed in the :ref:`sec_catalog` (or a fraction of a chromosome, using the
+``length_multiplier`` argument), keeping in mind that larger contigs will take
+longer to simulate. We could also specify a "generic" contig, which provides
+a segment of a given length with constant recombination rate, taken to be the
+average rate over all chromosomes for that species. Here, we define a contig
+of length 1 Mb:
+
+.. code-block:: python
+
+    contig = species.get_contig(length=1e6)
+    print(contig.recombination_map.sequence_length)
+    # 1000000.0
+    print(contig.recombination_map.mean_rate)
+    # 1.2313743222950562e-08
+    print(contig.mutation_rate)
+    # 1.29e-8
+
+The "sequence length" is the length in base pairs. Since we are using a generic
+contig, we cannot specify a recombination map so we get a "flat" map of
+constant recombination rate. The mutation rate defaults to the species average
+mutation rate, as no mutation rate was provided when defining the contig.
+
 Choose a sampling scheme, and simulate
 --------------------------------------
 
@@ -529,8 +572,8 @@ But, you can go crazy with the sample size!
     samples = model.get_samples(10)
     engine = stdpopsim.get_engine("msprime")
 
-Finally, we simulate the model with the contig length and number of samples we defined above.
-The simulation results are recorded in a tree sequence object
+Finally, we simulate the model with the contig length and number of samples we
+defined above. The simulation results are recorded in a tree sequence object
 (:class:`tskit.TreeSequence`).
 
 .. code-block:: python
@@ -551,17 +594,16 @@ Now, we do some simple checks that our simulation worked with
     print(ts.num_populations)
     # 1
     print(ts.num_mutations)
-    # 18187
+    # 1472
     print(ts.num_trees)
-    # 12203
+    # 1078
 
-As expected, there are 10 samples in one population. We can also see that it takes
-12203 distinct genealogical trees across these 5Mb of sequence,
-on which there were 18187 mutations
-(since we are not using a seed here, the number of mutations
-and trees will be slightly different for you). Try running the simulation again, and notice
-that the number of samples and populations stays the same, while the number of mutations
-and trees changes.
+As expected, there are 10 samples in one population. We can also see that it
+takes 1078 distinct genealogical trees across this 1Mb of sequence, on which
+there were 1472 mutations (since we are not using a seed here, the number of
+mutations and trees will be slightly different for you). Try running the
+simulation again, and notice that the number of samples and populations stays
+the same, while the number of mutations and trees changes.
 
 Output to VCF
 -------------
@@ -574,24 +616,24 @@ See the tskit documentation (:meth:`tskit.TreeSequence.write_vcf`) for more info
 .. code-block:: python
 
     with open("foo.vcf", "w") as vcf_file:
-        ts.write_vcf(vcf_file, contig_id="2")
+        ts.write_vcf(vcf_file, contig_id="0")
 
 Taking a look at the vcf file, we see something like this:
 
 .. code-block:: none
 
     ##fileformat=VCFv4.2
-    ##source=tskit 0.2.2
+    ##source=tskit 0.3.5
     ##FILTER=<ID=PASS,Description="All filters passed">
-    ##contig=<ID=2,length=12159969>
+    ##contig=<ID=0,length=1000000>
     ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
-    #CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	tsk_0	tsk_1	tsk_2	tsk_3	tsk_4	tsk_5	tsk_6	tsk_7	tsk_8	tsk_9
-    2	96	.	0	1	.	PASS	.	GT	0	0	1	0	1	0	0	0	1	0
-    2	129	.	0	1	.	PASS	.	GT	0	0	0	0	0	0	0	0	1	0
-    2	436	.	0	1	.	PASS	.	GT	0	0	0	0	0	1	0	0	0	0
-    2	466	.	0	1	.	PASS	.	GT	0	0	1	0	1	0	0	0	0	0
-    2	558	.	0	1	.	PASS	.	GT	0	0	0	0	0	0	0	0	1	0
-    2	992	.	0	1	.	PASS	.	GT	1	1	0	1	0	1	1	1	0	1
+    #CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  tsk_0   tsk_1   tsk_2   tsk_3   tsk_4   tsk_5   tsk_6   tsk_7   tsk_8   tsk_9
+    0       3608    .       A       T       .       PASS    .       GT      0       0       1       0       0       0       0       0       0       0
+    0       5598    .       T       G       .       PASS    .       GT      1       0       0       0       1       1       0       0       1       0
+    0       6190    .       C       A       .       PASS    .       GT      1       0       0       0       0       0       0       0       0       0
+    0       6479    .       C       T       .       PASS    .       GT      0       0       0       0       0       1       0       0       1       0
+    0       6556    .       A       C       .       PASS    .       GT      0       0       0       1       0       0       0       0       0       0
+    0       6648    .       T       A       .       PASS    .       GT      1       1       1       0       1       1       1       1       1       1
 
 
 Using the SLiM engine
@@ -635,9 +677,11 @@ it doesn't affect this very much either).
    import stdpopsim
 
    species = stdpopsim.get_species("HomSap")
-   contig = species.get_contig("chr22", length_multiplier=0.1)
-   # default is a flat genetic map
    model = species.get_demographic_model("Africa_1T12")
+   contig = species.get_contig(
+       "chr22", length_multiplier=0.1, mutation_rate=model.mutation_rate
+   )
+   # default is a flat genetic map with average rate across chr22
    samples = model.get_samples(200)
 
 
@@ -726,7 +770,9 @@ we'd do the following
 
 .. code-block:: python
 
-   contig = species.get_contig("chr22", genetic_map="HapMapII_GRCh37")
+   contig = species.get_contig(
+       "chr22", genetic_map="HapMapII_GRCh37", mutation_rate=model.mutation_rate
+   )
    samples = model.get_samples(200)
    engine = stdpopsim.get_engine("slim")
    ts = engine.simulate(model, contig, samples, slim_burn_in=0.1, slim_scaling_factor=10)
@@ -863,6 +909,7 @@ Using the ``--help-genetic-maps`` option, we can also see what recombination map
 are available:
 
 .. command-output:: stdpopsim HomSap --help-genetic-maps
+    :ellipsis: 20
 
 Let's go with ``HapMapII_GRCh37``.
 The next command simulates 4 samples of chromosome 1 from each of the four
@@ -873,7 +920,7 @@ To check that we have set up the simulation correctly, we may first wish to perf
 dry run using the ``-D`` option.
 This will print information about the simulation to the terminal:
 
-.. command-output:: stdpopsim HomSap -c chr1 -o   afr-america-chr1.trees -s 13 -g HapMapII_GRCh37 -d AmericanAdmixture_4B11 4 4 4 4 -D
+.. command-output:: stdpopsim HomSap -c chr1 -o afr-america-chr1.trees -s 13 -g HapMapII_GRCh37 -d AmericanAdmixture_4B11 4 4 4 4 -D
 
 Once we're sure, we can remove the ``-D`` flag to run the simulation. (Note: This took around 8 minutes to run on a laptop.)
 
diff --git a/stdpopsim/catalog/BosTau/demographic_models.py b/stdpopsim/catalog/BosTau/demographic_models.py
index 9129a027a..19e52f48e 100644
--- a/stdpopsim/catalog/BosTau/demographic_models.py
+++ b/stdpopsim/catalog/BosTau/demographic_models.py
@@ -23,7 +23,6 @@ def _HolsteinFriesian_1M13():
             reasons={stdpopsim.CiteReason.DEM_MODEL},
         )
     ]
-
     return stdpopsim.DemographicModel(
         id=id,
         description=description,
diff --git a/stdpopsim/catalog/HomSap/demographic_models.py b/stdpopsim/catalog/HomSap/demographic_models.py
index 4e5906853..0a31b9a20 100644
--- a/stdpopsim/catalog/HomSap/demographic_models.py
+++ b/stdpopsim/catalog/HomSap/demographic_models.py
@@ -137,7 +137,7 @@ 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
+        The model is derived from the Tennessen 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
@@ -292,7 +292,11 @@ def _america():
         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.
+        population of African ancestry, 1/3 European, and 1/2 Asian. Note that this
+        demographic model was not inferred, and the assumed mutation rate is smaller
+        than used in the inference from Gravel et al. Resulting levels of diversity
+        may be considerably lower than observed in human population genetic data for
+        these populations.
     """
     populations = [
         stdpopsim.Population(id="AFR", description="Contemporary African population"),
@@ -598,6 +602,7 @@ def _zigzag():
     ]
 
     generation_time = 30
+
     N0 = 5 * 14312
 
     g_1 = 0.02302578256
diff --git a/stdpopsim/cli.py b/stdpopsim/cli.py
index bf5fa0cda..a3838a10c 100644
--- a/stdpopsim/cli.py
+++ b/stdpopsim/cli.py
@@ -293,6 +293,11 @@ def write_citations(engine, model, contig, species):
     """
     cite_str = ["If you use this simulation in published work, please cite:"]
     for citation in get_citations(engine, model, contig, species):
+        if (
+            stdpopsim.citations.CiteReason.MUT_RATE in citation.reasons
+            and model.mutation_rate is not None
+        ):
+            continue
         cite_str.append("\n")
         cite_str.append(citation.displaystr())
     logger.warning("".join(cite_str))
@@ -518,6 +523,7 @@ def run_simulation(args):
             length=args.length,
             inclusion_mask=args.inclusion_mask,
             exclusion_mask=args.exclusion_mask,
+            mutation_rate=model.mutation_rate,
         )
         engine = stdpopsim.get_engine(args.engine)
         logger.info(
@@ -587,7 +593,17 @@ def write_simulation_summary(engine, model, contig, samples, seed=None):
     # Get information about relevant contig
     gmap = "None" if contig.genetic_map is None else contig.genetic_map.id
     mean_recomb_rate = contig.recombination_map.mean_rate
-    mut_rate = contig.mutation_rate
+    # use the model mutation rate, if provided with the demographic model
+    if model.mutation_rate is not None:
+        mut_rate = model.mutation_rate
+        logging.info(
+            f"using mutation rate from demographic model ({model.mutation_rate})"
+        )
+    else:
+        mut_rate = contig.mutation_rate
+        logging.info(
+            f"using mutation rate from species contig ({contig.mutation_rate})"
+        )
     contig_len = contig.recombination_map.sequence_length
     dry_run_text += "Contig Description:\n"
     dry_run_text += f"{indent}Contig length: {contig_len}\n"
diff --git a/stdpopsim/engines.py b/stdpopsim/engines.py
index 5f66d19be..4e229d7af 100644
--- a/stdpopsim/engines.py
+++ b/stdpopsim/engines.py
@@ -5,6 +5,7 @@
 import msprime
 import stdpopsim
 import numpy as np
+import math
 
 logger = logging.getLogger(__name__)
 
@@ -112,6 +113,19 @@ def _warn_zigzag(self, demographic_model):
                 "https://github.com/popsim-consortium/stdpopsim/issues/745"
             )
 
+    def _warn_mutation_rate_mismatch(self, contig, demographic_model):
+        if demographic_model.mutation_rate is not None and not math.isclose(
+            demographic_model.mutation_rate, contig.mutation_rate
+        ):
+            warnings.warn(
+                "The demographic model has mutation rate "
+                f"{demographic_model.mutation_rate}, but this simulation used the "
+                f"contig's mutation rate {contig.mutation_rate}. Diversity levels "
+                "may be different than expected for this species. For details see "
+                "documentation at "
+                "https://stdpopsim.readthedocs.io/en/latest/tutorial.html"
+            )
+
 
 class _MsprimeEngine(Engine):
     id = "msprime"  #:
@@ -226,6 +240,7 @@ def simulate(
 
         # TODO: remove this after a release or two. See #745.
         self._warn_zigzag(demographic_model)
+        self._warn_mutation_rate_mismatch(contig, demographic_model)
 
         rng = np.random.default_rng(seed)
         seeds = rng.integers(1, 2 ** 31 - 1, size=2)
diff --git a/stdpopsim/models.py b/stdpopsim/models.py
index a13110300..31aff4dd5 100644
--- a/stdpopsim/models.py
+++ b/stdpopsim/models.py
@@ -69,6 +69,9 @@ class DemographicModel:
     :ivar citations: A list of :class:`Citation`, that describe the primary
         reference(s) for the model.
     :vartype citations: list of :class:`Citation`
+    :ivar mutation_rate: The mutation rate associated with the demographic model.
+        If no mutation rate is given, the species' default mutation rate is used.
+    :vartype mutation_rate: float
     """
 
     def __init__(
@@ -86,11 +89,13 @@ def __init__(
         population_id_map=None,
         populations=None,
         model=None,
+        mutation_rate=None,
     ):
         self.id = id
         self.description = description
         self.long_description = long_description
         self.generation_time = 1 if generation_time is None else generation_time
+        self.mutation_rate = mutation_rate
         self.citations = [] if citations is None else citations
         self.qc_model = qc_model
 
@@ -172,7 +177,7 @@ def register_qc(self, qc_model):
             )
         if self.qc_model is not None:
             raise ValueError(f"QC model already registered for {self.id}.")
-        self.qc_model = qc_model.model
+        self.qc_model = qc_model
 
     def get_samples(self, *args):
         """
diff --git a/stdpopsim/qc/HomSap.py b/stdpopsim/qc/HomSap.py
index cfb834b7c..51c18a90f 100644
--- a/stdpopsim/qc/HomSap.py
+++ b/stdpopsim/qc/HomSap.py
@@ -247,7 +247,7 @@ def BrowningAmerica():
         id=id,
         description=id,
         long_description=id,
-        generation_time=_species.generation_time,
+        generation_time=25,
         populations=populations,
         # Set population sizes at T=0
         # pop0 is Africa, pop1 is Europe, pop2 is Asia, pop3 is admixed
diff --git a/stdpopsim/slim_engine.py b/stdpopsim/slim_engine.py
index 40668d92b..25e98f2de 100644
--- a/stdpopsim/slim_engine.py
+++ b/stdpopsim/slim_engine.py
@@ -1088,6 +1088,7 @@ def simulate(
 
         # TODO: remove this after a release or two. See #745.
         self._warn_zigzag(demographic_model)
+        self._warn_mutation_rate_mismatch(contig, demographic_model)
 
         run_slim = not slim_script
 
diff --git a/stdpopsim/species.py b/stdpopsim/species.py
index 97b2798f8..193f08e48 100644
--- a/stdpopsim/species.py
+++ b/stdpopsim/species.py
@@ -146,6 +146,7 @@ def get_contig(
         genetic_map=None,
         length_multiplier=1,
         length=None,
+        mutation_rate=None,
         inclusion_mask=None,
         exclusion_mask=None,
     ):
@@ -169,6 +170,8 @@ def get_contig(
             same chromosome-specific mutation and recombination rates.
             This option cannot currently be used in conjunction with the
             ``genetic_map`` argument.
+        :param float mutation_rate: The per-base mutation rate. If none is given,
+            the mutation rate defaults to the rate specified by species chromosomes.
         :param inclusion_mask: If specified, simulated genomes are subset to only
             inlude regions given by the mask. The mask can be specified by the
             path and file name of a bed file or as a list or array of intervals
@@ -211,10 +214,13 @@ def get_contig(
                     L_tot += chrom_data.length
                     r_tot += chrom_data.length * chrom_data.recombination_rate
                     u_tot += chrom_data.length * chrom_data.mutation_rate
-            u = u_tot / L_tot
+            if mutation_rate is None:
+                mutation_rate = u_tot / L_tot
             r = r_tot / L_tot
             recomb_map = msprime.RateMap.uniform(length, r)
-            ret = stdpopsim.Contig(recombination_map=recomb_map, mutation_rate=u)
+            ret = stdpopsim.Contig(
+                recombination_map=recomb_map, mutation_rate=mutation_rate
+            )
         else:
             if length is not None:
                 raise ValueError("Cannot specify sequence length for named contig")
@@ -255,9 +261,12 @@ def get_contig(
                 else:
                     exclusion_intervals = exclusion_mask
 
+            if mutation_rate is None:
+                mutation_rate = chrom.mutation_rate
+
             ret = stdpopsim.Contig(
                 recombination_map=recomb_map,
-                mutation_rate=chrom.mutation_rate,
+                mutation_rate=mutation_rate,
                 genetic_map=gm,
                 inclusion_mask=inclusion_intervals,
                 exclusion_mask=exclusion_intervals,


Test patch:
diff --git a/tests/test_models.py b/tests/test_models.py
index 268dd48e1..bd2d391f9 100644
--- a/tests/test_models.py
+++ b/tests/test_models.py
@@ -4,6 +4,7 @@
 import unittest
 import sys
 import textwrap
+import copy
 
 import numpy as np
 import msprime
@@ -30,8 +31,8 @@ def test_debug_runs(self):
 
     def test_simulation_runs(self):
         # With a recombination_map of None, we simulate a coalescent without
-        # recombination in msprime, with no mutation.
-        contig = stdpopsim.Contig()
+        # recombination in msprime, with mutation rate equal to rate from model.
+        contig = stdpopsim.Contig(mutation_rate=self.model.mutation_rate)
         # Generate vector with 2 samples for each pop with sampling enabled
         sample_count = []
         for p in self.model.populations:
@@ -62,10 +63,20 @@ class QcdCatalogDemographicModelTestMixin(CatalogDemographicModelTestMixin):
 
     def test_qc_model_equal(self):
         d1 = self.model.model
-        d2 = self.model.qc_model
+        d2 = self.model.qc_model.model
         d1.assert_equivalent(d2, rel_tol=1e-5)
         assert d1 != d2
 
+    def test_generation_time_match(self):
+        g1 = self.model.generation_time
+        g2 = self.model.qc_model.generation_time
+        self.assertEqual(g1, g2)
+
+    def test_mutation_rate_match(self):
+        u1 = self.model.mutation_rate
+        u2 = self.model.qc_model.mutation_rate
+        self.assertEqual(u1, u2)
+
 
 # Add model specific test classes, derived from one of the above.
 qc_test_classes = []
@@ -131,6 +142,8 @@ def test_all_instances(self, model):
         assert len(model.long_description) > 0
         assert len(model.citations) > 0
         assert model.generation_time > 0
+        if model.mutation_rate is not None:
+            assert model.mutation_rate > 0
         model.model.validate()
 
 
@@ -315,3 +328,34 @@ def test_zigzag_produces_warning(self):
         for engine in stdpopsim.all_engines():
             with self.assertWarnsRegex(UserWarning, "Zigzag_1S14"):
                 engine.simulate(model, contig, samples, dry_run=True)
+
+
+class TestMutationRates(unittest.TestCase):
+    def test_mutation_rate_warning(self):
+        species = stdpopsim.get_species("HomSap")
+        model = copy.deepcopy(species.get_demographic_model("OutOfAfrica_3G09"))
+        ###
+        # unneeded after mutation rates are included in catalog!!
+        model.mutation_rate = 2.35e-8
+        ###
+        contig = species.get_contig("chr22")
+        samples = model.get_samples(10, 10, 10)
+        for engine in stdpopsim.all_engines():
+            with self.assertWarnsRegex(UserWarning, "mutation"):
+                engine.simulate(model, contig, samples, dry_run=True)
+
+    def test_mutation_rate_match(self):
+        species = stdpopsim.get_species("HomSap")
+        model = copy.deepcopy(species.get_demographic_model("OutOfAfrica_3G09"))
+        ###
+        # unneeded after mutation rates are included in catalog!!
+        model.mutation_rate = 2.35e-8
+        ###
+        contig = species.get_contig("chr22")
+        self.assertNotEqual(model.mutation_rate, contig.mutation_rate)
+        contig = species.get_contig("chr22", mutation_rate=model.mutation_rate)
+        self.assertEqual(model.mutation_rate, contig.mutation_rate)
+        contig = species.get_contig(length=100)
+        self.assertNotEqual(model.mutation_rate, contig.mutation_rate)
+        contig = species.get_contig(length=100, mutation_rate=model.mutation_rate)
+        self.assertEqual(model.mutation_rate, contig.mutation_rate)