diff --git a/stdpopsim/qc/MusMus.py b/stdpopsim/qc/MusMus.py
index 10bbe85f0..91fa37d84 100644
--- a/stdpopsim/qc/MusMus.py
+++ b/stdpopsim/qc/MusMus.py
@@ -286,3 +286,42 @@ def QC_CastaneusIndia():
 _species.get_demographic_model("MusculusKorea_1F22").register_qc(QC_MusculusKorea())
 
 _species.get_demographic_model("CastaneusIndia_1F22").register_qc(QC_CastaneusIndia())
+
+
+def QC_GammaB21():
+    id = "QC-GammaB21"
+    # M. m. castaneus from https://www.biorxiv.org/content/10.1101/2021.06.10.447924v2
+    # 0-fold site model with parameters from Tables S1, S2
+    neutral = stdpopsim.MutationType()
+    gamma_shape = 0.18617976
+    # Using polyDFE, the DFE is estimated in terms of
+    # the scaled selection coefficient for deleterious mutations,
+    # 2Nesd, where sd is the reduction in fitness experienced by
+    # an individual homozygous for the mutation
+    gamma_mean = -50044.583
+    Ne = 420000  # effective population size from methods
+    gamma_mean /= Ne
+    h = 0.5
+    negative = stdpopsim.MutationType(
+        dominance_coeff=h,
+        distribution_type="g",
+        # PolyDFE gives two times the selection coefficient on a homozygote,
+        # so divide mean by two
+        distribution_args=[
+            np.around(gamma_mean / 2, decimals=4),
+            np.around(gamma_shape, decimals=3),
+        ],
+    )
+    # Tom's original code @ TBooker/MuridRodentProject/blob/master/bin/SlimFunctions.py
+    # parses the DFE assuming 1/3 of sites are neutral
+    prop_neutral = 0.334
+    return stdpopsim.DFE(
+        id=id,
+        description=id,
+        long_description=id,
+        mutation_types=[neutral, negative],
+        proportions=[prop_neutral, 1 - prop_neutral],
+    )
+
+
+_species.get_dfe("Gamma_B21").register_qc(QC_GammaB21())