From e536bfd9ce29faa31f30c448f3464af576f3f3d3 Mon Sep 17 00:00:00 2001
From: Ilan Gronau <ilan.gronau@idc.ac.il>
Date: Thu, 12 Oct 2023 19:45:29 +0300
Subject: [PATCH] Implemented 2-epoch demographic model for Vaquita

---
 .../PhoSin/Vaquita2Epoch_1R22.csv             |  5 ++
 stdpopsim/catalog/PhoSin/__init__.py          |  1 +
 .../catalog/PhoSin/demographic_models.py      | 57 +++++++++++++++++++
 3 files changed, 63 insertions(+)
 create mode 100644 docs/parameter_tables/PhoSin/Vaquita2Epoch_1R22.csv
 create mode 100644 stdpopsim/catalog/PhoSin/demographic_models.py

diff --git a/docs/parameter_tables/PhoSin/Vaquita2Epoch_1R22.csv b/docs/parameter_tables/PhoSin/Vaquita2Epoch_1R22.csv
new file mode 100644
index 000000000..0ceddcac3
--- /dev/null
+++ b/docs/parameter_tables/PhoSin/Vaquita2Epoch_1R22.csv
@@ -0,0 +1,5 @@
+Population size,"4,485",Ancestral pop. size
+Population size,"2,807",Pop. size during second epoch
+Epoch Time (gen.),"2,162",Start time of second epoch
+Generation time (yrs.),11.9,Average generation interval
+Mutation rate,5.83e-9,Per-base per-generation mutation rate
diff --git a/stdpopsim/catalog/PhoSin/__init__.py b/stdpopsim/catalog/PhoSin/__init__.py
index a518c0337..7e49c4342 100644
--- a/stdpopsim/catalog/PhoSin/__init__.py
+++ b/stdpopsim/catalog/PhoSin/__init__.py
@@ -2,3 +2,4 @@
 Catalog definitions for PhoSin (Ensembl ID='phocoena_sinus')
 """
 from . import species  # noqa: F401
+from . import demographic_models  # noqa: F401
diff --git a/stdpopsim/catalog/PhoSin/demographic_models.py b/stdpopsim/catalog/PhoSin/demographic_models.py
new file mode 100644
index 000000000..01568ffed
--- /dev/null
+++ b/stdpopsim/catalog/PhoSin/demographic_models.py
@@ -0,0 +1,57 @@
+import msprime
+import stdpopsim
+
+_species = stdpopsim.get_species("PhoSin")
+
+
+def _2epoch():
+    N_anc = 4485
+    N_curr = 2807
+    T = 2162
+    populations = [
+        stdpopsim.Population(
+            id="Vaquita",
+            description="Vaquita (Phocoena sinus)",
+        )
+    ]
+
+    return stdpopsim.DemographicModel(
+        id="Vaquita2Epoch_1R22",
+        description="Vaquita two epoch model",
+        long_description="""
+            A two-epoch demographic model estimated using dadi from site
+            frequency spectrum at putatively neutrally evolving regions of
+            the genome identified as those located >10 kb from coding
+            sequences which did not overlap with CpG islands.
+            Population genomic data obtained from 20 individuals sequenced
+            at mean coverage 60x.
+            Robinson et al. (2022) reports several inferred models in Supp
+            Table S2. This is the 2-epoch model inferred by dadi, which is
+            also depicted in Main Figure 1E.
+            Size changes from N_anc to N_curr in time T.
+        """,
+        populations=populations,
+        citations=[
+            stdpopsim.Citation(
+                author="Robinson et al.",
+                year=2022,
+                doi="https://doi.org/10.1126/science.abm1742",
+                reasons={stdpopsim.CiteReason.DEM_MODEL},
+            )
+        ],
+        generation_time=11.9,
+        mutation_rate=5.83e-9,
+        population_configurations=[
+            msprime.PopulationConfiguration(
+                initial_size=N_curr, metadata=populations[0].asdict()
+            )
+        ],
+        demographic_events=[
+            msprime.PopulationParametersChange(
+                time=T, initial_size=N_anc, population_id=0
+            )
+        ],
+    )
+
+
+_species.add_demographic_model(_2epoch())