Repository: tskit-dev/msprime (language: python)

Issue:
adding microsatellite mutation models 
i'm hoping to add some simple mutation models that would be appropriate for microsatellite data.
i think this will be quite useful to folks who are working on non-model organisms, ecological genetics, 
and maybe even forensics? 

The simplest, most classical model is the Stepwise Mutation Model of [Kimura and Ohta (1978)](https://dx.doi.org/10.1073%2Fpnas.75.6.2868). 
I reckon that's a good place to start, although [this nice paper](https://doi.org/10.1093/molbev/msg084) by Calabrese and Durrett 
lays out a number of other models that would be simple enough to implement. 

Here is how I would imagine the implementation to look for the stepwise microsat model. 
generally these models have a _lo_ and _hi_ parameter which governs the bounds on the
number of repeat units of the microsatelllite 

```python
class SMM_MicrosatMutationModel(MatrixMutationModel):
    """
    Stepwise mutation model (Kimura and Ohta, 1978), 
    for microsatellite repeats.

    This is a :class:`.MatrixMutationModel` with alleles ``[lo, .. , hi]``,
    root distribution is uniform ``(lo, hi)``, and a transition matrix that reflects
    single-step mutation.

    :param int lo: lower bound of the microsatellite repeat length
    :param int hi: upper bound of the microsatellite repeat length
        See :ref:`??` for more details.
    """

    def __init__(self, lo, hi):
        alleles = list(range(lo, hi+1))
        root_distribution = [1/len(alleles)] * len(alleles)
        transition_matrix = np.zeros((len(alleles), len(alleles)))
        for i in range(len(alleles)):
            if i == 0:
                transition_matrix[i, i+1] = 1
            elif i == len(alleles)-1:
                transition_matrix[i, i-1] = 1
            else:
                transition_matrix[i, i-1] = 0.5
                transition_matrix[i, i+1] = 0.5

        super().__init__(alleles, root_distribution, transition_matrix)
```

what do others think? 


Discussion / hints:
LGTM. I'd make ``hi`` exclusive just to keep consistent with things like ``range``, but otherwise looks very sensible.

So, to be concrete, the plan is to make alleles equal to ``[str(j) for j in range(lo, hi)]`` and then if you have allele j, you can only ever mutate to ``j - 1`` or ``j + 1``?
> So, to be concrete, the plan is to make alleles equal to `[str(j) for j in range(lo, hi)]` and then if you have allele j, you can only ever mutate to `j - 1` or `j + 1`?

yep that's the SMM model in a nutshell

This sounds fine, although we should make sure to tell people not to set `lo` and `hi` to the actual min and max observed in their data. Alternatively, we could modify the [infinite alleles](https://github.com/tskit-dev/msprime/blob/main/lib/mutgen.c#L415) or SLiM model code (not use the mutation matrix model), so that we can set `hi=Inf` (or, something large). It would be straightforward.
> Alternatively, we could modify the infinite alleles or SLiM model code (not use the mutation matrix model), so that we can set hi=Inf (or, something large). It would be straightforward.

That's not quite the same thing though is it? With infinite alleles you can only every go from j to j + 1, whereas the model above you can go back to j - 1 with equal proba.
> That's not quite the same thing though is it? With infinite alleles you can only every go from j to j + 1, whereas the model above you can go back to j - 1 with equal proba.

It's totally not the same thing - infinite alleles maintains a global counter for the next allele, and this would look at the previous state to choose a new one, but it provides a template. (so, maybe the SLiM mutation model code is better, as it actually modifies the previous allele)
Right, gotcha
> This sounds fine, although we should make sure to tell people not to set `lo` and `hi` to the actual min and max observed in their data. Alternatively, we could modify the [infinite alleles](https://github.com/tskit-dev/msprime/blob/main/lib/mutgen.c#L415) or SLiM model code (not use the mutation matrix model), so that we can set `hi=Inf` (or, something large). It would be straightforward.

another issue is the ancestral distribution. I think uniformly distributed is probably wrong?
The matrix method @andrewkern outlined above would work, though, right? Probably a good place to start if someone is doing this rather than diving into the C code immediately.
Yeah i think it should work just fine. 
Unless the number of alleles is 10_000 or something it may be good enough.
Good enough, but I don't think anyone's doing it right now, and it'd be nice to not clutter up this list of models with out of date ones?  Re: the root distribution - it'd be natural to allow a slight left/right bias and use the stationary distribution (I think there's a Slatkin paper on this?).
> Good enough, but I don't think anyone's doing it right now, and it'd be nice to not clutter up this list of models with out of date ones?

not sure what you mean? SMM is the most popular microsat model in analyses and simulations AFAIK
I mean if we implement another slightly different but better SMM later. But I suppose that allowing `hi=Inf` is backwards-compatible, as is adding left/right jump probs. So, never mind.

Fix patch:
diff --git a/.gitignore b/.gitignore
index ba70f9cc8..1e0d02841 100644
--- a/.gitignore
+++ b/.gitignore
@@ -93,3 +93,6 @@ tmp__NOBACKUP__
 
 # Mac stuff
 .DS_Store
+
+# vscode stuff
+.vscode
\ No newline at end of file
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7b222a13c..635f5712b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,15 @@
 # Changelog
 
+## [1.2.1] - 2022-XX-XX
+
+**New features**
+
+- Add a `MicrosatMutationModel` mutation model class, that
+  represents a generalized interface for constructing mutational
+  models appropriate to STRs. In addition 3 specific microsat models
+  are added `SMM`, `TPM`, and `EL2`.  
+  ({issue}`2013`, {user}`andrewkern`).
+
 ## [1.2.0] - 2022-05-18
 
 **New features**
@@ -10,7 +20,7 @@
 **Bug fixes**:
 
 - Fix rare assertion trip in the single sweep model caused by numerical jitter.
-  ({issue}`1966`, {pr}`2038`, {user}`jeromekelleher `, {user}`molpopgen`)
+  ({issue}`1966`, {pr}`2038`, {user}`jeromekelleher`, {user}`molpopgen`)
 
 - Fix edge case in `Demography.from_old_style()`
   ({issue}`2047`, {pr}`2048`, {user}`grahamgower`)
@@ -23,8 +33,9 @@
 ## [1.1.1] - 2022-02-10
 
 **Bug fixes**:
+
 - Fix (very) rare assertion trip caused by underlying GSL bug.
-  ({issue}`1997`, {pr}`2000`, {user}`chriscrsmith `, {user}`molpopgen`,
+  ({issue}`1997`, {pr}`2000`, {user}`chriscrsmith`, {user}`molpopgen`,
   {user}`andrewkern`)
 
 **Maintenance**:
@@ -94,7 +105,6 @@ This is a bugfix release recommended for all users.
 - Fix printing small values in rate maps ({pr}`1906`, {issue}`1905`,
   {user}`petrelharp`).
 
-
 ## [1.0.2] - 2021-06-29
 
 Improved Demes support and minor bugfixes.
@@ -126,7 +136,6 @@ Minor feature release with experimental Demes support.
 - Initial support for Demes via the ``Demography.from_demes`` method.
   ({pr}`1662`, {issue}`1675`, {user}`jeromekelleher`, {user}`apragsdale`)
 
-
 ## [1.0.0] - 2021-04-14
 
 Msprime 1.0 is a major update, recommended for all users. It introduces
@@ -331,6 +340,7 @@ Alpha release for testing the tskit/msprime split.
 Minor bugfix release.
 
 **New features**:
+
 - Add provenance recording option to simplify (#601)
 - Minor performance improvement (#598)
 
@@ -450,7 +460,7 @@ compatible, however.
 - The `SparseTree.mutations()` and `TreeSequence.mutations()` iterators no
   longer support tuple-like access to values. For example, code like
 
-  >   > for x, u, j in ts.mutations():
+  > > for x, u, j in ts.mutations():
   >
   > : print("mutation at position", x, "node = ", u)
   >
@@ -561,14 +571,14 @@ upcoming functionality.
 - The tuples returned by the `mutations` method contains an element.
   This will break code doing things like
 
-  >   > for pos, node in ts.mutations():
+  > > for pos, node in ts.mutations():
   >
   > : print(pos, node)
   >
   For better forward compatibility, code should use named attributes
   rather than positional access:
 
-  >   > for mutation in ts.mutations():
+  > > for mutation in ts.mutations():
   >
   > : print(mutation.position, mutation.node)
   >
@@ -794,4 +804,3 @@ populations.
 - Added a {}`get_pairwise_diversity` method to the TreeSequence
   class to efficiently calculate the population genetics
   statistic pi.
-
diff --git a/docs/api.md b/docs/api.md
index 85a0bde8b..abfa6bee8 100644
--- a/docs/api.md
+++ b/docs/api.md
@@ -12,8 +12,8 @@ for discussion and examples of individual features.
 ```{eval-rst}
 .. currentmodule:: msprime
 ```
-### Ancestry
 
+### Ancestry
 
 ```{eval-rst}
 .. autosummary::
@@ -44,6 +44,10 @@ for discussion and examples of individual features.
   PAM
   BinaryMutationModel
   MatrixMutationModel
+  MicrosatMutationModel
+  SMM
+  TPM
+  EL2
   InfiniteAlleles
   SLiMMutationModel
 ```
@@ -102,12 +106,10 @@ for discussion and examples of individual features.
 .. autoclass:: msprime.SmcPrimeApproxCoalescent
 ```
 
-
 ```{eval-rst}
 .. autoclass:: msprime.DiscreteTimeWrightFisher
 ```
 
-
 ```{eval-rst}
 .. autoclass:: msprime.FixedPedigree
 ```
@@ -126,7 +128,6 @@ for discussion and examples of individual features.
 
 ### Mutations
 
-
 ```{eval-rst}
 .. autofunction:: msprime.sim_mutations
 ```
@@ -141,6 +142,22 @@ for discussion and examples of individual features.
 .. autoclass:: msprime.MatrixMutationModel()
 ```
 
+```{eval-rst}
+.. autoclass:: msprime.MicrosatMutationModel()
+```
+
+```{eval-rst}
+.. autoclass:: msprime.SMM()
+```
+
+```{eval-rst}
+.. autoclass:: msprime.TPM()
+```
+
+```{eval-rst}
+.. autoclass:: msprime.EL2()
+```
+
 ```{eval-rst}
 .. autoclass:: msprime.BinaryMutationModel()
 ```
@@ -178,7 +195,6 @@ for discussion and examples of individual features.
 .. autoclass:: msprime.SLiMMutationModel()
 ```
 
-
 (sec_api_node_flags)=
 
 ### Node flags
@@ -227,7 +243,6 @@ The node was created by a census event. Please see the
 
 ### Rate maps
 
-
 ```{eval-rst}
 .. autoclass:: msprime.RateMap
     :members:
@@ -235,7 +250,6 @@ The node was created by a census event. Please see the
 
 ### Demography
 
-
 ```{eval-rst}
 .. autoclass:: msprime.Demography()
     :members:
diff --git a/docs/mutations.md b/docs/mutations.md
index cba3d7a40..a96a06c34 100644
--- a/docs/mutations.md
+++ b/docs/mutations.md
@@ -43,13 +43,16 @@ For example, here we simulate the ancestry of two diploids (and therefore
 four sample genomes) and store the resulting {class}`tskit.TreeSequence`
 object in the ``ts`` variable. We then simulate some mutations on
 that tree sequence, and show the resulting tree with mutations:
+
 ```{code-cell}
 ts = msprime.sim_ancestry(2, sequence_length=100, random_seed=1234)
 mts = msprime.sim_mutations(ts, rate=0.01, random_seed=5678)
 SVG(mts.draw_svg())
 ```
+
 We can see the variation data produced by these mutations most simply
 via the {meth}`tskit.TreeSequence.variants` method:
+
 ```{code-cell}
 for var in mts.variants():
     print(var.site.position, var.alleles, var.genotypes, sep="\t")
@@ -64,7 +67,6 @@ for more information.
 
 ## Quick reference
 
-
 {func}`.sim_mutations`
 : Add simulated mutations to a tree sequence
 
@@ -94,6 +96,15 @@ for more information.
 {class}`.MatrixMutationModel` (General finite state model)
 : Superclass of mutation models with a finite set of states
 
+{class}`.SMM` (Microsatellites)
+: Stepwise mutation model for microsatellites
+
+{class}`.TPM` (Microsatellites)
+: Equal rate, two-phase model for microsatellites
+
+{class}`.EL2` (Microsatellites)
+: Equal rate, two-phase, linear bias model for microsatellites
+
 {class}`.InfiniteAlleles` (Integers)
 : A generic infinite-alleles mutation model
 
@@ -153,8 +164,6 @@ position of the {class}`.RateMap` is the same as the
 
 :::
 
-
-
 (sec_mutations_discrete)=
 
 ## Discrete or continuous
@@ -181,7 +190,7 @@ mts.tables.sites.position
 :::{note}
 
 Using `discrete_genome=False` means that the mutation model will conform
-to the classic _infinite sites_ assumption, where each mutation in the simulation
+to the classic *infinite sites* assumption, where each mutation in the simulation
 occurs at a new site.
 
 :::
@@ -223,7 +232,6 @@ print(mts.tables.mutations.time)
 As explained in {ref}`the following section <sec_mutations_existing>`, reversing the order of these two
 lines will result in an error, as older mutations must be added first.
 
-
 (sec_mutations_silent)=
 
 ## Silent mutations
@@ -301,7 +309,6 @@ print(f"The tree sequence has {num_silent} silent mutations,\n"
 
 Indeed, 1182 - 75 = 1107.
 
-
 (sec_mutations_existing)=
 
 ## Existing mutations
@@ -372,7 +379,6 @@ with equal probability, which we can check by looking at the *transition matrix*
 whose rows give the probability that each allele mutates into each other allele
 (see {ref}`sec_mutations_matrix_mutations_models` for more):
 
-
 ```{code-cell}
 jc69 = msprime.JC69()
 jc69.transition_matrix
@@ -457,7 +463,6 @@ Here's the takeaways from this section:
     for more details.
 ```
 
-
 (sec_mutations_models)=
 
 ## Models
@@ -476,6 +481,9 @@ Here are the available models; they are documented in more detail below.
 - {class}`.BLOSUM62`: The BLOSUM62 model of time-reversible amino acid mutation
 - {class}`.PAM`: The PAM model of time-reversible amino acid mutation
 - {class}`.MatrixMutationModel`: Superclass of the specific mutation models with a finite set of states
+- {class}`.SMM`: Stepwise mutation model for microsatellite repeat copy number. Ohta and Kimura ('73)
+- {class}`.TPM`: Two-phase mutation model for microsatellite repeat copy number. DiRienzo et al. ('94)
+- {class}`.EL2`: Two-phase mutation model, equal rate, linear bias model for microsatellite repeat copy number. Garza et al. ('95)
 - {class}`.InfiniteAlleles`: A generic infinite-alleles mutation model
 - {class}`.SLiMMutationModel`: An infinite-alleles model of mutation producing SLiM-style mutations
 
@@ -611,7 +619,6 @@ Pt = scipy.linalg.expm(t * rate * Q)[i,j]
 If the top of a branch of length {math}`t` has allele {math}`i`,
 the bottom of the branch has allele {math}`j` with probability `Pt[i,j]`.
 
-
 (sec_mutations_adjusting_for_silent)=
 
 #### Adjusting mutation rates for silent mutations
@@ -661,9 +668,65 @@ print(f"Genetic diversity: {theta}.")
 
 That's pretty close to 0.001! The difference is within statistical error.
 
-
 (sec_mutations_mutation_infinite_alleles)=
 
+### Microsatellite Mutation Models
+
+We have implemented a number of mutational models that are commonly used
+to model the evolution of microsatellite repeats.
+The basic idea here is that the number of copies of a given repeat is
+tracked, and its evolution over time subject to one of a number of potential
+biases.
+
+Consider the [Combined DNA Index System](https://en.wikipedia.org/wiki/Combined_DNA_Index_System) (CODIS), a US National database
+of genotypes maintained by the justice system for use in forensic DNA analysis.
+The CODIS system relies on 20 core microsatellite loci, that are unlinked across
+the human genome. Let's simulate a sample of 20 individuals at 20 unlinked loci,
+and then apply the {class}`.SMM` model to the sample.
+
+Finally we will summarize variation across our loci, by
+iterating through the variants at each position / locus.
+We will print the position of the locus, then a list of allelic states, followed
+by haploid genotypes for each chromosome at that site.
+
+```{code-cell} python
+
+ts = msprime.sim_ancestry(20, random_seed=2, sequence_length=20,
+    recombination_rate=0.5, population_size=100_000)
+model = msprime.SMM(1, 10)
+mts = msprime.sim_mutations(ts, rate=1e-5, model=model, random_seed=1)
+for var in mts.variants():
+    print("Site at:", var.site.position)
+    print("    Alleles:", var.alleles)
+    print("    Genotypes:", var.genotypes)
+```
+
+We can also look at the dependence of the variance in repeat
+number on mutation rate under the SMM model. Much has been
+written over the years on how to estimate the mutation
+rate from the variance in repeat number
+(e.g., [Zhivotovsky and Feldman, 1995](https://pubmed.ncbi.nlm.nih.gov/8524801/)),
+so let's simulate that signal:
+
+```{code-cell} python
+from matplotlib import pyplot as plt
+
+rates = [1e-2,1e-3, 1e-4, 1e-5, 1e-6]
+model = msprime.SMM()
+vars = []
+ts = msprime.sim_ancestry(100, random_seed=2, sequence_length=1, population_size=100_000)
+for r in rates:
+    mts = msprime.sim_mutations(ts, rate=r, model=model, random_seed=1)
+    vars.append(mts.genotype_matrix().var())
+
+plt.plot(rates, vars)
+plt.xlabel('mutation rate')
+plt.ylabel('variance in repeat number')
+```
+
+So here we are seeing a clear relationship between the variance
+in repeat number and mutation rate, as expected. 🔥🔥🔥
+
 ### Infinite Alleles Mutation Models
 
 You can also use a model of *infinite alleles* mutation: where each new mutation produces a unique,
@@ -695,7 +758,6 @@ Apparently, there were 20 mutations at this site, but the alleles present in the
 Note that all other mutations on the tree are not observed in the population as they have
 been are "overwritten" by subsequent mutations.
 
-
 :::{warning}
 
 Neither this nor the next infinite alleles mutation model check to see if the alleles
@@ -776,7 +838,6 @@ SVG(t.draw_svg(mutation_labels=ml, node_labels={}, size=(400, 300)))
 
 ```
 
-
 Note what has happened here: on the top branch on the right side of the tree,
 with the first model we added two mutations: first a mutation with ID `0`,
 then a mutation with ID `3`.
@@ -789,4 +850,3 @@ This was already present in the tree sequence, so its derived state is not modif
 has "erased" the type 2 mutations `100` and `102`.
 If you want a different arrangement,
 you can go back and edit the derived states (and metadata) as you like.
-
diff --git a/msprime/__init__.py b/msprime/__init__.py
index eeeb60037..884ed24b4 100644
--- a/msprime/__init__.py
+++ b/msprime/__init__.py
@@ -70,7 +70,11 @@
     InfiniteAlleles,
     JC69,
     MatrixMutationModel,
+    MicrosatMutationModel,
     MutationModel,
+    SMM,
+    TPM,
+    EL2,
     NUCLEOTIDES,
     PAM,
     SLiMMutationModel,
@@ -123,6 +127,10 @@
     "JC69",
     "MassMigration",
     "MatrixMutationModel",
+    "MicrosatMutationModel",
+    "SMM",
+    "TPM",
+    "EL2",
     "MigrationRateChange",
     "MutationModel",
     "NODE_IS_CA_EVENT",
diff --git a/msprime/mutations.py b/msprime/mutations.py
index a195131fa..6090baab3 100644
--- a/msprime/mutations.py
+++ b/msprime/mutations.py
@@ -199,6 +199,315 @@ def __init__(self, state_independent=False):
         super().__init__(alleles, root_distribution, transition_matrix)
 
 
+def general_microsat_rate_matrix(s, u, v, p, m, lo, hi):
+    """
+    Construct the general microsat rate matrix associated
+    with the Sainudiin et al. (2004) parameterization.
+    While the original paper implements a continuous time markov chain,
+    here we implement the jump process rate matrix equivalent for consistency
+    with the :class:`.MatrixMutationModel`.
+
+    In particular if :math:`q_{ij}` is an element of the originally defined
+    transition matrix, we define :math:`M = \\max_i \\left( \\sum_j q_{ij} \\right)`.
+    Then, as a :class:`.MatrixMutationModel` this model has the transition
+    matrix :math:`P_{ij} = q_{ij} / M`, and :math:`P_{ii} = 1 - \\sum_{j\\neq i} P_{ij}`.
+
+    :param float s: strength of length dependence on mutation rate
+     (i.e. s=0 uniform rate)
+    :param float u: constant bias parameter; U ~ (0, 1)
+    :param float v: linear bias parameter; v ~ (-inf, inf)
+    :param float p: probability of a single step mutation
+    :param float m: success prob of truncated gamma distribution
+    :param int lo: lower bound on repeat number
+    :param int hi: upper bound on repeat number
+
+    :return: rate matrix
+
+    """
+
+    # Helper functions for the Sainudiin et al. rate matrix
+    def alpha_microsat(i, u, v, lo):
+        return max(0, min(1, u - (v * (i - lo))))
+
+    def beta_microsat(i, s, lo):
+        """
+        beta helper function from Sainudiin et al. (2004)
+        represents the mutation rate from allele i
+        :param float s: strength of length dependence on mutation rate
+        (i.e. s=0 uniform rate)
+        :param int lo: lower bound
+        :param int hi: upper bound
+        """
+        return 1 + ((i - lo) * s)
+
+    def gamma_microsat(m, i, j, lo, hi):
+        """
+        gamma helper function from Sainudiin et al. (2004)
+        :param float m: prob success
+        :param int i: in-state
+        :param int j: out-state
+        :param int lo: lower bound of microsat repeat length distribution
+        :param int hi: upper bound of microsat repeat length distribution
+        """
+        num = m * (1 - m) ** (np.abs(i - j) - 1)
+        if lo <= i < j <= hi:
+            denom = 1 - (1 - m) ** (hi - i)
+        elif lo <= j < i <= hi:
+            denom = 1 - (1 - m) ** (i - lo)
+        else:
+            raise TypeError("microsat model error. i and j must be between lo and hi")
+
+        return num / denom
+
+    nalleles = hi - lo + 1
+    rate_matrix = np.zeros((nalleles, nalleles))
+    for i in range(lo, hi + 1):
+        for j in range(lo, hi + 1):
+            # shift coords to account for 0-indexing in rate matrix
+            ii = i - lo
+            jj = j - lo
+            if i == j - 1:
+                rate_matrix[ii, jj] = (
+                    beta_microsat(i, s, lo)
+                    * alpha_microsat(i, u, v, lo)
+                    * (p + ((1 - p) * gamma_microsat(m, i, j, lo, hi)))
+                )
+            elif i < j - 1:
+                rate_matrix[ii, jj] = (
+                    beta_microsat(i, s, lo)
+                    * alpha_microsat(i, u, v, lo)
+                    * (1 - p)
+                    * gamma_microsat(m, i, j, lo, hi)
+                )
+            elif i == j + 1:
+                rate_matrix[ii, jj] = (
+                    beta_microsat(i, s, lo)
+                    * (1 - alpha_microsat(i, u, v, lo))
+                    * (p + ((1 - p) * gamma_microsat(m, i, j, lo, hi)))
+                )
+            elif i > j + 1:
+                rate_matrix[ii, jj] = (
+                    beta_microsat(i, s, lo)
+                    * (1 - alpha_microsat(i, u, v, lo))
+                    * (1 - p)
+                    * gamma_microsat(m, i, j, lo, hi)
+                )
+            else:
+                rate_matrix[ii, jj] = 0
+    # scale by max row sum; then normalize to 1
+    row_sums = rate_matrix.sum(axis=1, dtype="float64")
+    alpha = max(row_sums)
+    rate_matrix /= alpha
+    row_sums = rate_matrix.sum(axis=1, dtype="float64")
+    np.fill_diagonal(rate_matrix, 1.0 - row_sums)
+    return rate_matrix
+
+
+class MicrosatMutationModel(MatrixMutationModel):
+    """
+    General class for Microsatellite mutation models
+    which parameterizes the transition matrix describing
+    changes among alleles using 5 parameters according
+    to the model of Sainudiin et al. (2004).
+
+    Concretely, mutation rates under this model are computed as follows.
+    Let ``qij`` be the transition rate from equation (3) in Sainudiin et al,
+    and let ``M`` be the largest row sum of the matrix whose i,j-th entry
+    is ``qij`` and whose diagonals are zero. Then if the mutation rate
+    passed to :func:`.sim_mutations` is ``m``, the mutation rate from
+    allele ``i`` to allele `j`` is ``m * qij / M``.
+
+    :param float s: strength of length dependence on mutation rate
+        (i.e. s=0 uniform rate)
+    :param float u: constant bias parameter; must be between 0 and 1
+    :param float v: linear bias parameter
+    :param float p: probability of a single step mutation
+    :param float m: 1 / (mean multistep mutation size)
+    :param int lo: lower bound on repeat number
+    :param int hi: upper bound on repeat number
+    :param list[float] root_distribution: An array of float values of length hi-lo+1
+        which should sum to 1, and give the probabilities used to choose
+        the ancestral state of each microsatellite.
+        Defaults to the stationary distribution of the model.
+
+    """
+
+    def __init__(
+        self,
+        s=0,
+        u=0.5,
+        v=0.0,
+        p=1.0,
+        m=1.0,
+        lo=2,
+        hi=50,
+        root_distribution=None,
+    ):
+        if np.isinf(lo):
+            raise ValueError("lo is infinite")
+        if np.isinf(hi):
+            raise ValueError("hi is infinite")
+        if lo > hi:
+            raise ValueError("lo > hi")
+        if not lo >= 0:
+            raise ValueError("lo must be >= 0")
+        if s <= -1 / (hi - lo + 1):
+            raise ValueError("s must be > -1/(hi-lo+1)")
+        if u <= 0 or u >= 1:
+            raise ValueError("u must be between 0 and 1")
+        if type(lo) == float:
+            lo = int(lo)
+        if type(hi) == float:
+            hi = int(hi)
+        alleles = [str(int(x)) for x in range(lo, hi + 1)]
+        transition_matrix = general_microsat_rate_matrix(s, u, v, p, m, lo, hi)
+        if root_distribution is None:
+            # solve for stationary distribution
+            S, U = np.linalg.eig(transition_matrix.T)
+            U = np.real_if_close(U, tol=1)
+            stationary = np.array(U[:, np.where(np.abs(S - 1.0) < 1e-8)[0][0]])
+            stationary = stationary / np.sum(stationary)
+            root_distribution = stationary
+        super().__init__(alleles, root_distribution, transition_matrix)
+
+
+class SMM(MicrosatMutationModel):
+    """
+    Stepwise mutation model (Ohta and Kimura, 1978),
+    for microsatellite repeat number.
+
+    This is a :class:`.MicrosatMutationModel` with alleles ``[lo, .. , hi]``,
+    a root distribution (stationary distribution by default),
+    and a transition matrix that allows only mutations that
+    change the number of repeats by +/- 1: concretely, if the mutation
+    rate is ``m``, then an allele ``k`` mutates to ``k+1`` and ``k-1``
+    at rate ``m/2`` each,
+    except that mutations above ``hi`` or below ``lo`` result in no change.
+
+    Note that the values of ``lo`` and ``hi`` are in units of repeat number,
+    and are unrelated to repeat length
+    (e.g., the repeat itself can be dinucleotide, tri-nucleotide, etcetera).
+
+    :param int lo: lower bound of the microsatellite repeat copy number
+    :param int hi: upper bound of the microsatellite repeat copy number
+    :param list[float] root_distribution: An array of float values of length hi-lo+1
+        which should sum to 1, and give the probabilities used to choose
+        the ancestral state of each microsatellite.
+        Defaults to the value of stationary distribution
+        of the model.
+
+    """
+
+    def __init__(self, lo=2, hi=50, root_distribution=None):
+        # SMM specific parameters
+        s = 0
+        u = 0.5
+        v = 0
+        p = 0
+        m = 1
+        super().__init__(s, u, v, p, m, lo, hi, root_distribution)
+
+
+class TPM(MicrosatMutationModel):
+    """
+    Two-phase mutation model of DiRienzo et al., (1994),
+    https://doi.org/10.1073/pnas.91.8.3166
+
+    This models evolution of microsatellite repeat number in a manner
+    that allows for multi-step mutations in copy number of microsatellite repeats.
+    This is parameterized by `p` the probability of a single step mutation,
+    and `m` the success probability of the truncated gamma distribution
+    describing the distribution of longer steps such that the
+    mean multi-step mutation is length `1/m`.
+    For this model both `p` and `m` need to be set to < 1.0.
+
+    This is a :class:`.MicrosatMutationModel` with alleles ``[lo, .. , hi]``.
+    For a precise definition of the mutation rates, see :class:`.MicrosatMutationModel`
+    with ``s=0``, ``u=0.5``, and ``v=0``.
+    Default values of other parameters in the model
+    are those estimated in Sainudiin et al. (2004).
+
+    Note that the values of ``lo`` and ``hi`` are in units of repeat number,
+    and are unrelated to repeat length
+    (e.g., the repeat itself can be dinucleotide, tri-nucleotide, etcetera).
+
+    :param float p: probability of a single step mutation
+    :param float m: success prob of truncated gamma distribution
+    :param int lo: lower bound of the microsatellite repeat copy number
+    :param int hi: upper bound of the microsatellite repeat copy number
+    :param list[float] root_distribution: An array of float values of length hi-lo+1
+        which should sum to 1. These values are used to determine the ancestral state of
+        each microsatellite. Defaults to the stationary distribution
+        of the model
+
+    """
+
+    def __init__(self, p=0.9, m=0.93, lo=2, hi=50, root_distribution=None):
+        if not (0 < p < 1.0):
+            raise ValueError("p must be between 0 and 1")
+        if not (0 < m < 1.0):
+            raise ValueError("m must be between 0 and 1")
+        # TPM-specific fixed params
+        s = 0
+        u = 0.5
+        v = 0
+        super().__init__(s, u, v, p, m, lo, hi, root_distribution)
+
+
+class EL2(MicrosatMutationModel):
+    """
+    Equal rate, linear biased, two-phase mutation model of Garza et al., (1995),
+    [10.1093/oxfordjournals.molbev.a040239] as parameterized by
+    the Sainudiin et al. (2004) model.
+
+    This models evolution of microsatellite repeat number in a manner
+    that allows for multi-step mutations in copy number of microsatellite repeats.
+    This is parameterized by `p` the probability of a single step mutation,
+    and `m` the success probability of the truncated gamma distribution
+    describing the distribution of longer steps such that the
+    mean multi-step mutation is length `1/m`.
+    For this model both `p` and `m` need to be set to < 1.0.
+
+    In addition EL2 has a constant bias parameter `u`, where
+    :math:`0 < u < 1` and a linear bias parameter `v`, where
+    :math:`-\\infty < v < \\infty`.
+
+    This is a :class:`.MicrosatMutationModel` with alleles ``[lo, .. , hi]``
+    For a precise definition of the mutation rates, see :class:`.MicrosatMutationModel`
+    with ``s=0``.
+    Default values for other parameters of the model
+    are those estimated in Sainudiin et al. (2004).
+
+    Note that the values of ``lo`` and ``hi`` are in units of repeat number,
+    and are unrelated to repeat length
+    (e.g., the repeat itself can be dinucleotide, tri-nucleotide, etcetera).
+
+    :param float p: probability of a single step mutation
+    :param float m: success prob of truncated gamma distribution
+    :param float u: constant bias parameter
+    :param float v: linear bias parameter
+    :param int lo: lower bound of the microsatellite repeat copy number
+    :param int hi: upper bound of the microsatellite repeat copy number
+    :param list[float] root_distribution: An array of float values of length hi-lo+1
+        which should sum to 1. These values are used to determine the ancestral state of
+        each microsatellite. Defaults to the stationary distribution
+        of the model
+
+    """
+
+    def __init__(self, m=0.43, u=0.68, v=0.037, lo=2, hi=50, root_distribution=None):
+        if not (0 < m < 1.0):
+            raise ValueError("m must be between 0 and 1")
+        if not (0 < u < 1.0):
+            raise ValueError("u must be between 0 and 1")
+        if not (-np.inf < v < np.inf):
+            raise ValueError("v must be between -inf and inf")
+        # EL2 specific params
+        p = 0.0
+        s = 0
+        super().__init__(s, u, v, p, m, lo, hi, root_distribution)
+
+
 class JC69(MatrixMutationModel):
     """
     Jukes-Cantor mutation model (Jukes and Cantor 1969). Based on the standard ACGT


Test patch:
diff --git a/tests/test_mutations.py b/tests/test_mutations.py
index 412f562d6..ab8ddfb47 100644
--- a/tests/test_mutations.py
+++ b/tests/test_mutations.py
@@ -339,6 +339,150 @@ def test_new_model(self):
         self.validate_model(model)
         self.validate_stationary(model)
 
+    def test_SMM(self):
+        model = msprime.SMM(lo=2, hi=50)
+        # test for symmetry of SMM trans mat
+        for i in range(model.transition_matrix.shape[0]):
+            for j in range(i):
+                if i != j:
+                    assert (
+                        model.transition_matrix[i, j] == model.transition_matrix[j, i]
+                    )
+                if np.abs(i - j) > 1:
+                    assert model.transition_matrix[i, j] == 0.0
+                elif np.abs(i - j) == 1:
+                    assert (
+                        model.transition_matrix[i, j] == model.transition_matrix[0, 1]
+                    )
+
+        self.validate_model(model)
+        self.validate_stationary(model)
+
+    def test_TPM(self):
+        lo = 2
+        hi = 50
+        m = 0.9
+        p = 0.9
+        model = msprime.TPM(p=p, m=m, lo=lo, hi=hi)
+        self.validate_model(model)
+        self.validate_stationary(model)
+
+        # test for ratio of TPM trans mat vals
+        for i in range(lo, hi + 1):
+            if lo < i < hi:
+                ii = i - lo
+                exp = (p + (m / (1 - (1 - m) ** (hi - i))) * (1 - p)) / (
+                    p + (m / (1 - (1 - m) ** (i - lo)) * (1 - p))
+                )
+                obs = (
+                    model.transition_matrix[ii, ii + 1]
+                    / model.transition_matrix[ii, ii - 1]
+                )
+                assert np.isclose(exp, obs)
+
+    def test_EL2(self):
+        model = msprime.EL2()
+        self.validate_model(model)
+        self.validate_stationary(model)
+
+        lo = 1
+        hi = 10
+        m = 0.43
+        u = 0.68
+        v = 0.037
+        p = 0
+        model = msprime.EL2(lo=lo, hi=hi)
+        for i in range(lo, hi + 1):
+            if lo < i < hi:
+                ii = i - lo
+                g_ratio = (p + (m / (1 - (1 - m) ** (hi - i))) * (1 - p)) / (
+                    p + (m / (1 - (1 - m) ** (i - lo)) * (1 - p))
+                )
+                a_ratio = u - (v * (i - lo))
+                a_ratio /= 1 - a_ratio
+                exp = g_ratio * a_ratio
+                obs = (
+                    model.transition_matrix[ii, ii + 1]
+                    / model.transition_matrix[ii, ii - 1]
+                )
+                assert np.abs(exp - obs) < 1e-6
+                # summation tests
+                sum_ij1 = model.transition_matrix[ii, :ii].sum()
+                sum_ij2 = model.transition_matrix[ii, ii:].sum()
+                exp = sum_ij2 / sum_ij1
+                assert np.abs(exp - a_ratio) < 1e-6
+
+    @pytest.mark.parametrize(
+        "m, u, v, p, s",
+        [
+            (1.0, 0.62, 0.015, 0, 0.88),
+            (0.43, 0.68, 0.037, 0, 0.88),
+            (0.43, 0.68, 0.037, 0.5, 0.88),
+            (0.43, 0.68, 0.037, 1.0, 0.88),
+        ],
+    )
+    def test_MicrosatMutationModel(self, m, u, v, p, s):
+        lo = 1
+        hi = 10
+        model = msprime.MicrosatMutationModel(s, u, v, p, m, lo, hi)
+        row_sums = model.transition_matrix.sum(axis=1)
+        d = row_sums - model.transition_matrix.diagonal()
+        betas = [1 + ((i - lo) * s) for i in range(lo, hi + 1)]
+        alphas = [u - (v * (i - lo)) for i in range(lo, hi + 1)]
+        for j in range(len(betas)):
+            exp = betas[0] / betas[j] * alphas[0] / (1 - alphas[j])
+            assert ((d[0] / d[j]) - exp) < 1e-6
+        for j in range(2, len(betas) - 1):
+            exp = betas[1] / betas[j]
+            assert ((d[0] / d[j]) - exp) < 1e-6
+
+    @pytest.mark.parametrize(
+        "m, u, v, p, s, lo, hi",
+        [
+            (1.0, 0.62, 0.015, 0, 0.88, 1, np.inf),
+            (0.43, 0.68, 0.037, 0, 0.88, np.inf, 10),
+            (0.43, 0.68, 0.037, 0, 0.88, 10, 1),
+            (0.43, 0.68, 0.037, 0, 0.88, -1, 4),
+            (0.43, 0.68, 0.037, 1.0, -10, 1, 4),
+            (0.43, 10, 0.037, 0, 0.88, 1, 5),
+        ],
+    )
+    def test_bad_MicrosatMutationModel(self, m, u, v, p, s, lo, hi):
+        with pytest.raises(ValueError):
+            msprime.MicrosatMutationModel(s, u, v, p, m, lo, hi)
+
+    @pytest.mark.parametrize(
+        "m, u, v, p, s, lo, hi",
+        [
+            (1.0, 0.62, 0.015, 0, 0.88, 1, 10.0),
+            (0.43, 0.68, 0.037, 1.0, 0.5, 3.0, 10),
+            (0.83, 0.58, 0.37, 0, 0.88, 1, 20),
+            (0.03, 0.18, 0.5, 1.0, 1, 1, 4),
+        ],
+    )
+    def test_good_params_MicrosatMutationModel(self, m, u, v, p, s, lo, hi):
+        model = msprime.MicrosatMutationModel(s, u, v, p, m, lo, hi)
+        self.validate_model(model)
+        self.validate_stationary(model)
+
+    def test_bad_TPM_params(self):
+        with pytest.raises(ValueError):
+            msprime.TPM(1.1, 0.9, 2, 9)
+        with pytest.raises(ValueError):
+            msprime.TPM(0.9, 1.1, 2, 9)
+
+    def test_bad_EL2_params(self):
+        with pytest.raises(ValueError):
+            msprime.EL2(1.1, 0.9, 0.1, 0.1, 2, 9)
+        with pytest.raises(ValueError):
+            msprime.EL2(0.9, 1.1, 0.1, 0.1, 2, 9)
+        with pytest.raises(ValueError):
+            msprime.EL2(0.9, 0.9, -0.1, 0.1, 2, 9)
+        with pytest.raises(ValueError):
+            msprime.EL2(0.9, 0.9, 0.1, np.inf, 2, 9)
+        with pytest.raises(ValueError):
+            msprime.EL2(0.9, 0.9, np.inf, np.inf, 2, 9)
+
 
 class TestMutateRateMap:
     def test_no_muts(self):