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

Issue:
merge the msprime-1.0 feature branch into main
Should be fairly straight forward to do this once #775 goes in.
Merge branch 'msprime-1.0' (genomic element types) into main
Closes #776.


Discussion / hints:
@mufernando - I think the main changes in that branch are yours. Do you think we can make a PR to bring the changes into main now?
I'm planning to have look at merging this later today.
thanks, @grahamgower!
# [Codecov](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/877?src=pr&el=h1) Report
> Merging [#877](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/877?src=pr&el=desc) (2945027) into [main](https://codecov.io/gh/popsim-consortium/stdpopsim/commit/0f6824e556f3b6e1e6450b4a01546cc0c9396fc4?el=desc) (0f6824e) will **decrease** coverage by `0.01%`.
> The diff coverage is `99.23%`.

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

```diff
@@            Coverage Diff             @@
##             main     #877      +/-   ##
==========================================
- Coverage   99.38%   99.37%   -0.02%     
==========================================
  Files          62       62              
  Lines        2450     2551     +101     
  Branches      285      312      +27     
==========================================
+ Hits         2435     2535     +100     
  Misses          7        7              
- Partials        8        9       +1     
```


| [Impacted Files](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/877?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [stdpopsim/ext/selection.py](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/877/diff?src=pr&el=tree#diff-c3RkcG9wc2ltL2V4dC9zZWxlY3Rpb24ucHk=) | `100.00% <ø> (ø)` | |
| [stdpopsim/genomes.py](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/877/diff?src=pr&el=tree#diff-c3RkcG9wc2ltL2dlbm9tZXMucHk=) | `99.29% <98.79%> (-0.71%)` | :arrow_down: |
| [stdpopsim/slim\_engine.py](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/877/diff?src=pr&el=tree#diff-c3RkcG9wc2ltL3NsaW1fZW5naW5lLnB5) | `98.81% <100.00%> (+0.02%)` | :arrow_up: |
| [stdpopsim/utils.py](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/877/diff?src=pr&el=tree#diff-c3RkcG9wc2ltL3V0aWxzLnB5) | `100.00% <100.00%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/877?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/877?src=pr&el=footer). Last update [0f6824e...2945027](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/877?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

> @grahamgower, I need your input in some of these comments. Also, I guess I should go ahead and fix a lot of the stuff I pointed out. How do I do that in your PR?

Ok, cool. I should probably have just left this for you to do. You can pull down this PR following [the tskit dev docs](https://tskit.dev/tskit/docs/stable/development.html#checking-out-someone-else-s-pull-request). You *might* be able to push back to my branch (from which the PR was made), and your changes will just appear here. I'm not sure how this all works with the permissions though... If you get permission denied errors, just open a new PR and we can close this one.
Oh, and just a heads up that there might be conflicts with #839 that need to be resolved (here or there, depending on who wins the merge race). Ping @apragsdale.
> Oh, and just a heads up that there might be conflicts with #839 that need to be resolved (here or there, depending on who wins the merge race). Ping @apragsdale.

Might be biased here, but the changes to the species class are pretty minimal in #839. Since it's a bug fix I'd like to get it merged as quick as possible, and also so I can also get to updating mutation rates in the catalog itself. I don't see any _major_ conflicts with the genomic elements introduced here, so hopefully it would be straightforward to resolve here?
> the changes to the species class are pretty minimal in #839.

Ah, I see now that the changes to the Contig class have been removed from there (or I imagined that). Carry on...
ok, I'll finish this PR. but probably not until next week. I'm busy with #ProbGen21 stuff this week.
Hey @mufernando, this PR came up in the PopSim call today. Any chance you could take another look? Feel free to rebase/remerge/whatever, reopen it as a new pr, or whatever you think is easiest!
Yes, I can merge this PR tomorrow. Sorry for the delay, this slipped my
mind.
-- 
*Murillo Fernando Rodrigues*


On Tue, May 4, 2021 at 9:27 AM Graham Gower ***@***.***>
wrote:

> Hey @mufernando <https://github.com/mufernando>, this PR came up in the
> PopSim call today. Any chance you could take another look? Feel free to
> rebase/remerge/whatever, reopen it as a new pr, or whatever you think is
> easiest!
>
> —
> You are receiving this because you were mentioned.
> Reply to this email directly, view it on GitHub
> <https://github.com/popsim-consortium/stdpopsim/pull/877#issuecomment-832075012>,
> or unsubscribe
> <https://github.com/notifications/unsubscribe-auth/ACL5HFDY2GUTBYJ4E2DCRITTMAOALANCNFSM42UQSQXQ>
> .
>

I'll do this in #930.

Fix patch:
diff --git a/docs/selection_example.py b/docs/selection_example.py
index 5705b4c39..11e308022 100644
--- a/docs/selection_example.py
+++ b/docs/selection_example.py
@@ -1,6 +1,6 @@
 import random
 import logging
-
+import numpy as np
 import stdpopsim
 
 logger = logging.getLogger(__name__)
@@ -27,7 +27,7 @@ def adaptive_introgression(seed):
     # SLiM phase of the simulation has completed.
     positive = stdpopsim.ext.MutationType(convert_to_substitution=False)
     mutation_types = [positive]
-    mut_id = len(mutation_types)
+    mut_id = len(mutation_types) - 1
 
     # We need some demographic model parameters to set bounds on the timing
     # of random variables and extended_events (below).
@@ -149,7 +149,6 @@ def adaptive_introgression(seed):
         contig,
         samples,
         seed=rng.randrange(1, 2 ** 32),
-        mutation_types=mutation_types,
         extended_events=extended_events,
         slim_scaling_factor=10,
         slim_burn_in=0.1,
@@ -165,17 +164,17 @@ def KimDFE():
     Return neutral and negative MutationType()s representing a human DFE.
     Kim et al. (2018), p.23, http://doi.org/10.1371/journal.pgen.1007741
     """
-    neutral = stdpopsim.ext.MutationType(weight=1.0)
+    neutral = stdpopsim.ext.MutationType()
     gamma_shape = 0.186  # shape
     gamma_mean = -0.01314833  # expected value
     h = 0.5 / (1 - 7071.07 * gamma_mean)  # dominance coefficient
     negative = stdpopsim.ext.MutationType(
-        weight=2.31,
         dominance_coeff=h,
         distribution_type="g",
         distribution_args=[gamma_mean, gamma_shape],
     )
-    return [neutral, negative]
+    # note neutral mutations have 0 proportion because they are not simulated by SLiM
+    return {"mutation_types": [neutral, negative], "proportions": [0.0, 0.7]}
 
 
 def OutOfAfrica_3G09_with_DFE(seed):
@@ -187,7 +186,41 @@ def OutOfAfrica_3G09_with_DFE(seed):
     contig = species.get_contig("chr1", length_multiplier=0.001)
     samples = model.get_samples(100, 100, 100)  # YRI, CEU, CHB
 
-    mutation_types = KimDFE()
+    # neutral and deleterious mutations occur across the whole contig
+    contig.add_genomic_element_type(
+        intervals=np.array([[0, int(contig.length)]]), **KimDFE()
+    )
+
+    # Simulate.
+    engine = stdpopsim.get_engine("slim")
+    ts = engine.simulate(
+        model,
+        contig,
+        samples,
+        seed=seed,
+        slim_scaling_factor=10,
+        slim_burn_in=10,
+        # Set slim_script=True to print the script instead of running it.
+        # slim_script=True,
+    )
+    return ts
+
+
+def gene_with_noncoding_OutOfAfrica_3G09(seed):
+    """
+    Simulating a 1kb gene flanked by 1kb neutral regions. Within genes,
+    30% of the total influx of mutations are neutral and 70% are deleterious,
+    with the DFE from Kim et al. The HomSap/OutOfAfrica_3G09 model was simulated.
+    """
+    species = stdpopsim.get_species("HomSap")
+    model = species.get_demographic_model("OutOfAfrica_3G09")
+    contig = species.get_contig(length=3000)
+    samples = model.get_samples(100, 100, 100)  # YRI, CEU, CHB
+
+    # within the gene, KimDFE is used, outside genomic elements
+    # neutral muts are added with msprime
+    gene_interval = np.array([[1000, 2000]])
+    contig.add_genomic_element_type(intervals=gene_interval, **KimDFE())
 
     # Simulate.
     engine = stdpopsim.get_engine("slim")
@@ -196,7 +229,6 @@ def OutOfAfrica_3G09_with_DFE(seed):
         contig,
         samples,
         seed=seed,
-        mutation_types=mutation_types,
         slim_scaling_factor=10,
         slim_burn_in=10,
         # Set slim_script=True to print the script instead of running it.
@@ -222,3 +254,5 @@ def OutOfAfrica_3G09_with_DFE(seed):
     ts, T_mut, T_sel, s = adaptive_introgression(seed)
 
     ts = OutOfAfrica_3G09_with_DFE(seed)
+
+    ts = gene_with_noncoding_OutOfAfrica_3G09(seed)
diff --git a/stdpopsim/ext/selection.py b/stdpopsim/ext/selection.py
index c9e5e7c9f..96abb355d 100644
--- a/stdpopsim/ext/selection.py
+++ b/stdpopsim/ext/selection.py
@@ -13,15 +13,8 @@ class MutationType(object):
     distribution_type = attr.ib(default="f", type=str)
     distribution_args = attr.ib(factory=lambda: [0], type=list)
     convert_to_substitution = attr.ib(default=True, type=bool)
-    # MutationTypes with non-zero weight will be simulated by SLiM,
-    # using rates obtained from the relative weights of the types.
-    # I.e. the weights will be used in the ``proportion`` parameter
-    # to SLiM's :func:`initializeGenomicElementType()`.
-    weight = attr.ib(default=0, type=float)
 
     def __attrs_post_init__(self):
-        if self.weight < 0:
-            raise ValueError("weight must be >= 0.")
         if self.dominance_coeff < 0:
             raise ValueError(f"Invalid dominance coefficient {self.dominance_coeff}.")
 
@@ -63,19 +56,6 @@ def __attrs_post_init__(self):
         }[self.distribution_type]
 
 
-def slim_mutation_frac(mutation_types):
-    """
-    The fraction of mutations that should be added by SLiM.
-    The remainder are added by msprime.mutate() once the SLiM part
-    of the simulation is complete.
-    """
-    if mutation_types is None:
-        weighted = False
-    else:
-        weighted = any(mut_type.weight > 0 for mut_type in mutation_types)
-    return 1 if weighted else 0
-
-
 @attr.s
 class GenerationAfter(object):
     """
diff --git a/stdpopsim/genomes.py b/stdpopsim/genomes.py
index 85a710e6a..d1e93caa8 100644
--- a/stdpopsim/genomes.py
+++ b/stdpopsim/genomes.py
@@ -1,7 +1,13 @@
 """
-Infrastructure for defining information about species' genomes.
+Infrastructure for defining information about species' genomes and genomic
+regions to be simulated.
 """
+import logging
 import attr
+import numpy as np
+import stdpopsim
+
+logger = logging.getLogger(__name__)
 
 
 @attr.s
@@ -133,12 +139,41 @@ class Chromosome:
     synonyms = attr.ib(factory=list, kw_only=True)
 
 
+@attr.s(kw_only=True)
+class GenomicElementType(object):
+    """
+    Class that represents a map between genomic intervals and a group of
+    mutation types and the rate in which they occur within those intervals.
+
+    :ivar mutation_type_ids: Indices for the mutation types
+        (:class:`.ext.MutationType`) that can occur in genomic elements of this
+        type.
+    :vartype mutation_type_ids: list
+    :ivar proportions: Proportions (out of the mutation rate) in which each
+        mutation type can happen in genomic intervals of this type.
+    :vartype: proportions: numpy.array, dtype=np.float32
+    :ivar intervals: (n, 2) numpy array with [left, right) genomic intervals
+        that will be of this GenomicElementType. The intervals are used to
+        specify SLiM GenomicElements. Intervals must be non-overlapping.
+    :vartype intervals: numpy.array, dtype=np.int32
+    """
+
+    mutation_type_ids = attr.ib(factory=lambda: [], type=list)
+    proportions = attr.ib(factory=lambda: np.array([]), type=np.ndarray)
+    intervals = attr.ib(default=None, type=np.array)
+
+    def __attrs_post_init__(self):
+        _check_mut_proportions(self.proportions, self.mutation_type_ids)
+        stdpopsim.utils.check_intervals_validity(self.intervals)
+
+
 @attr.s
 class Contig:
     """
     Class representing a contiguous region of genome that is to be
     simulated. This contains the information about mutation rates
     and recombination rates that are needed to simulate this region.
+    Genomic element types can be provided.
 
     :ivar mutation_rate: The rate of mutation per base per generation.
     :vartype mutation_rate: float
@@ -153,6 +188,13 @@ class Contig:
     :ivar exclude: If True, ``mask_intervals`` specify regions to exclude. If False,
         ``mask_intervals`` specify regions in keep.
     :vartype exclude: bool
+    :ivar genomic_element_types: a list of :class:`.GenomicElementType` objects.
+        By default, the only element is a single genomic element type that spans
+        the entire contiguous region.
+    :vartype genomic_element_types: list
+    :ivar mutation_types: a list of :class:`.ext.MutationType` object. By default,
+        the only mutation type is neutral. See :class:`.ext.MutationType` for
+        more details.
 
     .. note::
         To run stdpopsim simulations with alternative, user-specified mutation
@@ -174,16 +216,117 @@ class Contig:
     genetic_map = attr.ib(default=None, kw_only=True)
     inclusion_mask = attr.ib(default=None, kw_only=True)
     exclusion_mask = attr.ib(default=None, kw_only=True)
+    genomic_element_types = attr.ib(factory=list, type=list, kw_only=True)
+    mutation_types = attr.ib(factory=list, type=list, kw_only=True)
+
+    @property
+    def all_intervals_array(self):
+        """
+        Returns an (n, 3) numpy array with the intervals across all genomic
+        element types in the form (left, right, type).
+        """
+        all_intervals = np.concatenate(
+            [
+                np.column_stack(
+                    (g.intervals[:, :2], np.full((g.intervals.shape[0]), i))
+                )
+                for i, g in enumerate(self.genomic_element_types)
+            ],
+            axis=0,
+        )
+        all_intervals = stdpopsim.utils.build_intervals_array(all_intervals)
+        return all_intervals
+
+    def add_mutation_types(self, mutation_types, proportions, genomic_element_type_id):
+        """
+        Adds mutation types with their respective proportions to the genomic
+        element type provided.
+        """
+        if genomic_element_type_id >= len(self.genomic_element_types):
+            raise ValueError(
+                f"Genomic element type of id {genomic_element_type_id} does not exist."
+            )
+
+        def _is_index(y, X):
+            for i, e in enumerate(X):
+                if e is y:
+                    return i
+            return -1
+
+        for mut_type, prop in zip(mutation_types, proportions):
+            mid = _is_index(mut_type, self.mutation_types)
+            if mid == -1:
+                self.mutation_types.append(mut_type)
+                mid = len(self.mutation_types) - 1
+            self.genomic_element_types[
+                genomic_element_type_id
+            ].mutation_type_ids.append(mid)
+            self.genomic_element_types[genomic_element_type_id].proportions = np.append(
+                self.genomic_element_types[genomic_element_type_id].proportions, [prop]
+            )
+        genomic_element_type = self.genomic_element_types[genomic_element_type_id]
+        _check_mut_proportions(
+            genomic_element_type.proportions, genomic_element_type.mutation_type_ids
+        )
+
+    def add_genomic_element_type(self, intervals, mutation_types, proportions):
+        if self.recombination_map is None:
+            raise ValueError(
+                "You cannot add genomic element types to a "
+                "contig that has no recombination map."
+            )
+        stdpopsim.utils.check_intervals_validity(intervals)
+        ge_type = stdpopsim.GenomicElementType(intervals=intervals)
+        self.genomic_element_types.append(ge_type)
+        getid = len(self.genomic_element_types) - 1
+        self.add_mutation_types(
+            mutation_types, proportions, genomic_element_type_id=getid
+        )
+
+    def clear_genomic_mutation_types(self):
+        self.genomic_element_types = []
+        self.mutation_types = []
+
+    def fully_neutral(self, slim_mutations=False, convert_to_substitution=True):
+        self.clear_genomic_mutation_types()
+        if self.recombination_map is None:
+            raise ValueError("You must specify a recombination map.")
+        self.add_genomic_element_type(
+            intervals=np.array([[0, int(self.recombination_map.sequence_length)]]),
+            mutation_types=[
+                stdpopsim.ext.MutationType(
+                    convert_to_substitution=convert_to_substitution
+                )
+            ],
+            proportions=[1 if slim_mutations else 0],
+        )
 
     def __str__(self):
         gmap = "None" if self.genetic_map is None else self.genetic_map.id
         s = (
             "Contig(length={:.2G}, recombination_rate={:.2G}, "
-            "mutation_rate={:.2G}, genetic_map={})"
+            "mutation_rate={:.2G}, genetic_map={}, genomic_element_types={}, "
+            "mutation_types={})"
         ).format(
             self.recombination_map.sequence_length,
             self.recombination_map.mean_rate,
             self.mutation_rate,
             gmap,
+            self.genomic_element_types,
+            self.mutation_types,
         )
         return s
+
+
+def _check_mut_proportions(proportions, mutation_type=None, atol=1e-15):
+    if mutation_type is not None and len(proportions) != len(mutation_type):
+        raise ValueError(
+            "The list of proportions and mutation types should be of the same length."
+        )
+    if any([not (0 <= prop <= 1) for prop in proportions]):
+        raise ValueError("Proportions must lie within the [0,1] interval.")
+    if sum(proportions) > 1 + atol:
+        raise ValueError(
+            "The sum of the proportions within genomic element types cannot"
+            " exceed 1."
+        )
diff --git a/stdpopsim/slim_engine.py b/stdpopsim/slim_engine.py
index 25e98f2de..d40b9f7af 100644
--- a/stdpopsim/slim_engine.py
+++ b/stdpopsim/slim_engine.py
@@ -56,6 +56,7 @@
 import numpy as np
 import msprime
 import pyslim
+import tskit
 
 logger = logging.getLogger(__name__)
 
@@ -73,8 +74,6 @@
 
     defineConstant("burn_in", $burn_in);
     defineConstant("generation_time", $generation_time);
-    defineConstant("mutation_rate", Q * $mutation_rate);
-    defineConstant("chromosome_length", $chromosome_length);
     defineConstant("trees_file", "$trees_file");
     defineConstant("pop_names", $pop_names);
 
@@ -88,8 +87,6 @@
     defineConstant("N", asInteger(_N/Q));
 
     initializeTreeSeq();
-    initializeMutationRate(mutation_rate);
-    initializeGenomicElement(g1, 0, chromosome_length-1);
     initializeRecombinationRate(recombination_rates, recombination_ends);
 }
 
@@ -521,18 +518,72 @@ def msprime_rm_to_slim_rm(recombination_map):
     return rates, ends[1:]
 
 
+def get_slim_fractions(contig):
+    return np.array([sum(g.proportions) for g in contig.genomic_element_types])
+
+
+def compute_msp_mutation_rate_map(contig, slim_fractions):
+    """
+    breaks = [0,10,20] means [0,10) and [10,20)
+    """
+    length = contig.recombination_map.sequence_length
+    breaks = [0]
+    rates = []
+    for (start, end, t) in contig.all_intervals_array:
+        if start not in breaks:
+            breaks.append(start)
+            rates.append(contig.mutation_rate)
+        breaks.append(end)
+        rates.append(
+            np.around((1 - slim_fractions[t]) * contig.mutation_rate, decimals=16)
+        )
+    if not np.isclose(breaks[-1], int(length)):
+        breaks.append(int(length))
+        rates.append(contig.mutation_rate)
+    assert len(breaks) == len(rates) + 1
+    return msprime.RateMap(position=breaks, rate=rates)
+
+
+def compute_slim_mutation_breaks_and_rates(contig, msp_mutation_rate_map):
+    """
+    breaks = [10,20] means [0,10] and (10,20]
+    """
+    # making sure the mut types are still concordant with the ge types
+    breaks = msp_mutation_rate_map.position
+    rates = msp_mutation_rate_map.rate
+    assert breaks[0] == 0
+    slim_breaks = [b - 1 for b in breaks[1:]]
+    slim_rates = [contig.mutation_rate - r for r in rates]
+    return (slim_breaks, slim_rates)
+
+
+def get_msp_and_slim_mutation_rate_maps(contig):
+    """
+    Returns a :class:`msprime.RateMap` with the mutation rate map for overlay
+    of neutral mutations and a tuple with the breakpoints and rates of mutations
+    for SLiM.
+    """
+    slim_fractions = get_slim_fractions(contig)
+    msp_mutation_rate_map = compute_msp_mutation_rate_map(contig, slim_fractions)
+
+    slim_breaks, slim_rates = compute_slim_mutation_breaks_and_rates(
+        contig, msp_mutation_rate_map
+    )
+    return (msp_mutation_rate_map, (slim_breaks, slim_rates))
+
+
 def slim_makescript(
     script_file,
     trees_file,
     demographic_model,
     contig,
     samples,
-    mutation_types,
     extended_events,
     scaling_factor,
     burn_in,
+    slim_rate_map,
 ):
-
+    mutation_types = contig.mutation_types
     pop_names = [pop.name for pop in demographic_model.model.populations]
     # Use copies of these so that the time frobbing below doesn't have
     # side-effects in the caller's model.
@@ -645,14 +696,10 @@ def fix_time(event):
         if hasattr(ee, "mutation_type_id"):
             mt_id = getattr(ee, "mutation_type_id")
             cls_name = ee.__class__.__name__
-            if mutation_types is None:
-                raise ValueError(
-                    f"Invalid {cls_name} event. No mutation types defined."
-                )
-            if not (0 < ee.mutation_type_id <= len(mutation_types)):
-                # FIXME: use zero-based indexes
+            if not (0 <= mt_id < len(mutation_types)):
                 raise ValueError(
-                    f"Invalid {cls_name} event with mutation type id {mt_id}."
+                    f"Invalid {cls_name} {ee.mutation_type_id} {len(mutation_types)} "
+                    f"event with mutation type id {mt_id}.{contig}"
                 )
         if hasattr(ee, "start_time") and hasattr(ee, "end_time"):
             # Now that GenerationAfter times have been accounted for, we can
@@ -662,6 +709,7 @@ def fix_time(event):
             stdpopsim.ext.validate_time_range(start_time, end_time)
 
         if isinstance(ee, stdpopsim.ext.DrawMutation):
+            # TODO: warning if mutation type id points to a neutral mut.
             time = ee.time * demographic_model.generation_time
             save = 1 if ee.save else 0
             drawn_mutations.append(
@@ -708,7 +756,7 @@ def fix_time(event):
                 cls_name = ee.__class__.__name__
                 raise ValueError(
                     f"Invalid {cls_name} event. No drawn mutation for "
-                    "mutation type id {ee.mutation_type_id}"
+                    f"mutation type id {ee.mutation_type_id}"
                 )
 
     printsc = functools.partial(print, file=script_file)
@@ -757,10 +805,8 @@ def fix_time(event):
         string.Template(_slim_upper).substitute(
             scaling_factor=scaling_factor,
             burn_in=float(burn_in),
-            chromosome_length=int(contig.recombination_map.sequence_length),
             recombination_rates=recomb_rates_str,
             recombination_ends=recomb_ends_str,
-            mutation_rate=contig.mutation_rate,
             generation_time=demographic_model.generation_time,
             trees_file=trees_file,
             pop_names=f"c({pop_names_str})",
@@ -802,29 +848,37 @@ def matrix2str(
 
         return "".join(s)
 
-    if mutation_types is None:
-        mutation_types = [stdpopsim.ext.MutationType()]
-
-    # Mutation type; genomic elements.
-    # FIXME: Change this to use zero-based indices---one-based indices are
-    #        inconsistent with everything else, e.g. population IDs.
-    for i, m in enumerate(mutation_types, 1):
+    # Mutation types.
+    for i, m in enumerate(mutation_types):
         distrib_args = [str(arg) for arg in m.distribution_args]
         distrib_args[m.Q_scaled_index] = "Q * " + distrib_args[m.Q_scaled_index]
         distrib_args = ", ".join(distrib_args)
         printsc(
-            f'    initializeMutationType("m{i}", {m.dominance_coeff}, '
+            f"    initializeMutationType({i}, {m.dominance_coeff}, "
             + f'"{m.distribution_type}", {distrib_args});'
         )
         if not m.convert_to_substitution:
             # T is the default for WF simulations.
             printsc(f"    m{i}.convertToSubstitution = F;")
-    mut_weights = ", ".join(str(m.weight) for m in mutation_types)
-    printsc(
-        '    initializeGenomicElementType("g1", '
-        + f"seq(1, {len(mutation_types)}), c({mut_weights}));"
-    )
-    printsc()
+    # Genomic element types.
+    genomic_element_types = contig.genomic_element_types
+    for j, g in enumerate(genomic_element_types):
+        mut_props = ", ".join([str(prop) for prop in g.proportions])
+        mut_types = ", ".join([str(mid) for mid in g.mutation_type_ids])
+        element_starts = np.array2string(g.intervals[:, 0], separator=", ")[1:-1]
+        # stdpopsim intervals are 0-based left inclusive, right exclusive, but
+        # SLiM intervals are right inclusive
+        element_ends = np.array2string((g.intervals[:, 1] - 1), separator=", ")[1:-1]
+        printsc(
+            f"    initializeGenomicElementType({j}, c({mut_types}), c({mut_props}));"
+        )
+        printsc(
+            f"    initializeGenomicElement({j}, c({element_starts}), c({element_ends}));"
+        )
+    # Mutation rate map.
+    mut_breaks = ", ".join([str(int(b)) for b in slim_rate_map[0]])
+    mut_rates = ", ".join([str(r) for r in slim_rate_map[1]])
+    printsc(f"    initializeMutationRate(c({mut_rates})*Q, c({mut_breaks}));")
 
     # Epoch times.
     printsc("    // Time of epoch boundaries, in years before present.")
@@ -1034,7 +1088,6 @@ def simulate(
         samples,
         *,
         seed=None,
-        mutation_types=None,
         extended_events=None,
         slim_path=None,
         slim_script=False,
@@ -1085,6 +1138,18 @@ def simulate(
                     "compare results across different values of the scaling factor."
                 )
             )
+        if contig.recombination_map is None:
+            raise ValueError(
+                "Your Contig object is not compatible with the SLiM "
+                "engine. You must specify a recombination map."
+            )
+        # if no GET were defined, fully neutral sim with overlay of muts using msprime
+        if len(contig.genomic_element_types) == 0:
+            contig.fully_neutral()
+
+        # figuring out mutations in SLiM and neutral mutations to be added
+        # by msprime later
+        msp_rate_map, slim_rate_map = get_msp_and_slim_mutation_rate_maps(contig)
 
         # TODO: remove this after a release or two. See #745.
         self._warn_zigzag(demographic_model)
@@ -1092,17 +1157,6 @@ def simulate(
 
         run_slim = not slim_script
 
-        # Ensure only "weighted" mutations are introduced by SLiM.
-        mutation_rate = contig.mutation_rate
-        slim_frac = stdpopsim.ext.slim_mutation_frac(mutation_types)
-        contig = stdpopsim.Contig(
-            recombination_map=contig.recombination_map,
-            mutation_rate=slim_frac * mutation_rate,
-            genetic_map=contig.genetic_map,
-            inclusion_mask=contig.inclusion_mask,
-            exclusion_mask=contig.exclusion_mask,
-        )
-
         mktemp = functools.partial(tempfile.NamedTemporaryFile, mode="w")
 
         @contextlib.contextmanager
@@ -1121,10 +1175,10 @@ def script_file_f():
                 demographic_model,
                 contig,
                 samples,
-                mutation_types,
                 extended_events,
                 slim_scaling_factor,
                 slim_burn_in,
+                slim_rate_map,
             )
 
             script_file.flush()
@@ -1139,10 +1193,10 @@ def script_file_f():
             if dry_run:
                 return None
 
-            ts = pyslim.load(ts_file.name)
+            ts = tskit.load(ts_file.name)
 
         ts = self._recap_and_rescale(
-            ts, seed, recap_epoch, contig, mutation_rate, slim_frac, slim_scaling_factor
+            ts, seed, recap_epoch, contig, slim_scaling_factor, msp_rate_map
         )
 
         if contig.inclusion_mask is not None:
@@ -1218,14 +1272,7 @@ def _simplify_remembered(self, ts):
         return ts.simplify(samples=list(nodes), filter_populations=False)
 
     def _recap_and_rescale(
-        self,
-        ts,
-        seed,
-        recap_epoch,
-        contig,
-        mutation_rate,
-        slim_frac,
-        slim_scaling_factor,
+        self, ts, seed, recap_epoch, contig, slim_scaling_factor, msp_rate_map
     ):
         """
         Apply post-SLiM transformations to ``ts``. This rescales node times,
@@ -1239,7 +1286,8 @@ def _recap_and_rescale(
         ts_metadata = tables.metadata
         ts_metadata["SLiM"]["generation"] *= slim_scaling_factor
         tables.metadata = ts_metadata
-        ts = pyslim.SlimTreeSequence.load_tables(tables)
+        ts = tables.tree_sequence()
+        slim_gen = ts.metadata["SLiM"]["generation"]
 
         rng = random.Random(seed)
         s1, s2 = rng.randrange(1, 2 ** 32), rng.randrange(1, 2 ** 32)
@@ -1250,38 +1298,52 @@ def _recap_and_rescale(
             )
             for pop in recap_epoch.populations
         ]
-        ts = ts.recapitate(
-            recombination_rate=contig.recombination_map.mean_rate,
+
+        demography = msprime.Demography.from_old_style(
             population_configurations=population_configurations,
             migration_matrix=recap_epoch.migration_matrix,
+        )
+
+        # TODO: using recombination map in recapitate. maps should be discrete.
+        ts = msprime.sim_ancestry(
+            initial_state=ts,
+            demography=demography,
+            recombination_rate=contig.recombination_map.mean_rate,
             random_seed=s1,
         )
 
         ts = self._simplify_remembered(ts)
 
-        if slim_frac < 1:
-            # Add mutations to SLiM part of trees.
-            rate = (1 - slim_frac) * mutation_rate
-            ts = pyslim.SlimTreeSequence(
-                msprime.mutate(
-                    ts,
-                    rate=rate,
-                    keep=True,
-                    random_seed=s2,
-                    end_time=ts.slim_generation,
-                )
-            )
+        # Figuring out SLiM mutation id metadata
+        max_id = -1
+        for mut in ts.mutations():
+            for d in mut.derived_state.split(","):
+                max_id = max(max_id, int(d))
+
+        # Creating a SLiM mutation model
+        model = msprime.SLiMMutationModel(
+            type=len(contig.mutation_types), next_id=max_id + 1
+        )
 
         # Add mutations to recapitated part of trees.
+        ts = msprime.sim_mutations(
+            ts,
+            rate=contig.mutation_rate,
+            model=model,
+            start_time=slim_gen,
+            keep=True,
+            random_seed=s2,
+        )
+
+        # Adding neutral mutations to the SLiM part
         s3 = rng.randrange(1, 2 ** 32)
-        ts = pyslim.SlimTreeSequence(
-            msprime.mutate(
-                ts,
-                rate=mutation_rate,
-                keep=True,
-                random_seed=s3,
-                start_time=ts.slim_generation,
-            )
+        ts = msprime.sim_mutations(
+            ts,
+            rate=msp_rate_map,
+            model=model,
+            end_time=slim_gen,
+            keep=True,
+            random_seed=s3,
         )
 
         return ts
@@ -1292,7 +1354,6 @@ def recap_and_rescale(
         demographic_model,
         contig,
         samples,
-        mutation_types=None,
         extended_events=None,
         slim_scaling_factor=1.0,
         seed=None,
@@ -1312,7 +1373,7 @@ def recap_and_rescale(
         generate the SLiM script with :func:`simulate`.
 
         :param ts: The tree sequence output by SLiM.
-        :type ts: :class:`pyslim.SlimTreeSequence`
+        :type ts: :class:`tskit.TreeSequence`
 
         .. warning::
             The :func:`recap_and_rescale` function is provided in the hope that
@@ -1321,14 +1382,15 @@ def recap_and_rescale(
             should be consulted to determine if the behaviour is appropriate
             for your case.
         """
-        # Only "weighted" mutations are introduced by SLiM.
-        mutation_rate = contig.mutation_rate
-        slim_frac = stdpopsim.ext.slim_mutation_frac(mutation_types)
-        contig = stdpopsim.Contig(
-            recombination_map=contig.recombination_map,
-            mutation_rate=slim_frac * mutation_rate,
-            genetic_map=contig.genetic_map,
-        )
+        if contig.recombination_map is None:
+            raise ValueError(
+                "Your Contig object is not compatible with the SLiM "
+                "engine. You must specify a recombination map."
+            )
+        # if no GET were defined, fully neutral sim with overlay of muts using msprime
+        if len(contig.genomic_element_types) == 0:
+            contig.fully_neutral()
+        msp_rate_map, slim_rate_map = get_msp_and_slim_mutation_rate_maps(contig)
 
         with open(os.devnull, "w") as script_file:
             recap_epoch = slim_makescript(
@@ -1337,14 +1399,14 @@ def recap_and_rescale(
                 demographic_model,
                 contig,
                 samples,
-                mutation_types,
                 extended_events,
                 slim_scaling_factor,
                 1,
+                slim_rate_map,
             )
 
         ts = self._recap_and_rescale(
-            ts, seed, recap_epoch, contig, mutation_rate, slim_frac, slim_scaling_factor
+            ts, seed, recap_epoch, contig, slim_scaling_factor, msp_rate_map
         )
         return ts
 
diff --git a/stdpopsim/utils.py b/stdpopsim/utils.py
index bf1756ef1..0bf3a6745 100644
--- a/stdpopsim/utils.py
+++ b/stdpopsim/utils.py
@@ -9,6 +9,7 @@
 import tarfile
 import contextlib
 import numpy as np
+import tskit
 
 
 def is_valid_demographic_model_id(model_id):
@@ -170,3 +171,35 @@ def add_if_unique(chrom, syn):
         # Mt
         if chrom.id == "MT":
             add_if_unique(chrom, "chr" + "M")
+
+
+def check_intervals_array_shape(intervals):
+    if len(intervals.shape) != 2 or intervals.shape[1] < 2:
+        raise ValueError(
+            "Intervals must be 2D objects with at least 2 columns " "[left, right)."
+        )
+
+
+def check_intervals_validity(intervals, start=0, end=np.inf):
+    check_intervals_array_shape(intervals)
+    if np.any(intervals[:, 0] >= intervals[:, 1]):
+        raise ValueError(
+            "Left positions cannot be greater than or equal to right positions."
+        )
+    if np.any(intervals[:, :2] > end) or np.any(intervals[:, :2] < start):
+        raise ValueError(f"Intervals must be within [{start}, {end}]")
+    if not np.all(intervals[1:, 0] >= intervals[:-1, 1]):
+        raise ValueError("Intervals must be non-overlapping.")
+
+
+def build_intervals_array(intervals, start=0, end=np.inf):
+    """
+    Converts a 2D list or numpy.array to an np.int32 array, which is sorted by
+    the first axis. It also checks for the validity of intervals and non-overlappingness.
+    """
+    intervals = tskit.util.safe_np_int_cast(intervals, dtype=np.int64)
+    check_intervals_array_shape(intervals)
+    sorter = intervals[:, 0].argsort()
+    intervals = intervals[sorter]
+    check_intervals_validity(intervals, start, end)
+    return intervals


Test patch:
diff --git a/tests/test_genomes.py b/tests/test_genomes.py
new file mode 100644
index 000000000..50fbcfc3f
--- /dev/null
+++ b/tests/test_genomes.py
@@ -0,0 +1,89 @@
+"""
+Tests for classes that hold information about genomic region to be simulated.
+"""
+import unittest
+import stdpopsim
+import numpy as np
+
+
+class TestGenomicElementType(unittest.TestCase):
+    def test_bad_proportions(self):
+        # no proportion but mut type, prop > 1, sum(prop) > 1
+        proportions = ([], [2], [0.5, 0.5, 0.5])
+        mut_type_ids = ([0], [0], [0])
+        for props in proportions:
+            with self.assertRaises(ValueError):
+                stdpopsim.GenomicElementType(
+                    intervals=np.array([[0, 1]]),
+                    mutation_type_ids=mut_type_ids,
+                    proportions=props,
+                )
+
+
+class TestContig(unittest.TestCase):
+    def get_test_contig(self):
+        species = stdpopsim.get_species("HomSap")
+        contig = species.get_contig("chr22", length_multiplier=0.001)
+        return contig
+
+    def test_all_intervals_array(self):
+        contig = self.get_test_contig()
+        contig.fully_neutral()
+        self.assertTrue(len(contig.all_intervals_array) == 1)
+
+        truth = np.array([[20, 30, 1], [30, 50, 0]])
+        contig.clear_genomic_mutation_types()
+        contig.add_genomic_element_type(
+            intervals=np.array([[30, 50, 0]]),
+            mutation_types=[stdpopsim.ext.MutationType()],
+            proportions=np.array([0]),
+        )
+        contig.add_genomic_element_type(
+            intervals=np.array([[20, 30, 1]]),
+            mutation_types=[stdpopsim.ext.MutationType()],
+            proportions=np.array([0]),
+        )
+        self.assertTrue((contig.all_intervals_array == truth).all())
+
+    def test_fully_neutral_error(self):
+        contig = self.get_test_contig()
+        contig.recombination_map = None
+        with self.assertRaises(ValueError):
+            contig.fully_neutral()
+
+    def test_add_mutation_type_errors(self):
+        contig = self.get_test_contig()
+        # diff length mutation types and proportions
+        with self.assertRaises(ValueError):
+            contig.add_mutation_types([stdpopsim.ext.MutationType()], [], 0)
+        # invalid ge type id
+        with self.assertRaises(ValueError):
+            contig.add_mutation_types([stdpopsim.ext.MutationType()], [0], 0)
+        contig.fully_neutral()
+        # sum props > 1
+        with self.assertRaises(ValueError):
+            contig.add_mutation_types([stdpopsim.ext.MutationType()], [1.1], 0)
+
+    def test_add_genomic_element_type_errors(self):
+        contig = self.get_test_contig()
+        contig.recombination_map = None
+        # can't add ge type without rec map
+        with self.assertRaises(ValueError):
+            contig.add_genomic_element_type(np.array([]), [], [])
+        contig = self.get_test_contig()
+        # bad intervals
+        with self.assertRaises(ValueError):
+            contig.add_genomic_element_type(np.array([10, 20]), [], [])
+
+    def test_add_genomic_element_type(self):
+        contig = self.get_test_contig()
+        self.assertTrue(contig.mutation_types == [])
+        self.assertTrue(contig.genomic_element_types == [])
+        intervals1 = np.array([[10, 20], [20, 40], [100, 200]])
+        intervals2 = np.array([[5, 10], [45, 50], [200, 250]])
+        m1 = stdpopsim.ext.MutationType()
+        m2 = stdpopsim.ext.MutationType()
+        contig.add_genomic_element_type(intervals1, [m1, m2], [0.5, 0.3])
+        contig.add_genomic_element_type(intervals2, [m1], [0.2])
+        self.assertTrue(len(contig.genomic_element_types) == 2)
+        self.assertTrue(len(contig.mutation_types) == 2)
diff --git a/tests/test_slim_engine.py b/tests/test_slim_engine.py
index 310761b81..84f3002c3 100644
--- a/tests/test_slim_engine.py
+++ b/tests/test_slim_engine.py
@@ -10,6 +10,7 @@
 import tempfile
 import math
 from unittest import mock
+import numpy as np
 
 import tskit
 import pyslim
@@ -114,6 +115,20 @@ def test_script_generation(self):
         )
         self.assertTrue("sim.registerLateEvent" in out)
 
+    def test_no_recombination_map(self):
+        engine = stdpopsim.get_engine("slim")
+        species = stdpopsim.get_species("HomSap")
+        contig = species.get_contig("chr1", genetic_map="HapMapII_GRCh37")
+        model = stdpopsim.PiecewiseConstantSize(species.population_size)
+        contig.recombination_map = None
+        with self.assertRaises(ValueError):
+            engine.simulate(
+                demographic_model=model,
+                contig=contig,
+                samples=model.get_samples(10),
+                dry_run=True,
+            )
+
     def test_recombination_map(self):
         engine = stdpopsim.get_engine("slim")
         species = stdpopsim.get_species("HomSap")
@@ -147,22 +162,19 @@ def test_simulate(self):
     def test_recap_and_rescale(self):
         engine = stdpopsim.get_engine("slim")
         species = stdpopsim.get_species("HomSap")
-        contig = species.get_contig("chr22", length_multiplier=0.001)
         model = species.get_demographic_model("OutOfAfrica_3G09")
         samples = model.get_samples(10, 10, 10)
-
-        for weight, seed in zip((0, 1), (1234, 2345)):
-            if weight:
-                mutation_types = [stdpopsim.ext.MutationType(weight=1)]
+        for proportion, seed in zip((0, 1), (1234, 2345)):
+            contig = species.get_contig("chr22", length_multiplier=0.001)
+            if proportion:
+                contig.fully_neutral(slim_mutations=True)
                 extended_events = None
             else:
-                mutation_types = None
                 extended_events = []
             ts1 = engine.simulate(
                 demographic_model=model,
                 contig=contig,
                 samples=samples,
-                mutation_types=mutation_types,
                 extended_events=extended_events,
                 slim_scaling_factor=10,
                 slim_burn_in=0,
@@ -172,7 +184,6 @@ def test_recap_and_rescale(self):
                 demographic_model=model,
                 contig=contig,
                 samples=samples,
-                mutation_types=mutation_types,
                 extended_events=extended_events,
                 slim_scaling_factor=10,
                 slim_burn_in=0,
@@ -183,7 +194,6 @@ def test_recap_and_rescale(self):
                 demographic_model=model,
                 contig=contig,
                 samples=samples,
-                mutation_types=mutation_types,
                 extended_events=extended_events,
                 slim_scaling_factor=10,
                 seed=seed,
@@ -518,6 +528,20 @@ def test_warning_when_scaling(self):
             with self.assertWarns(stdpopsim.SLiMScalingFactorWarning):
                 capture_output(stdpopsim.cli.stdpopsim_main, cmd.split())
 
+    def test_slim_engine_without_recombination_map(self):
+        engine, species, contig = self.triplet()
+        contig.recombination_map = None
+        model = stdpopsim.PiecewiseConstantSize(1000)
+        samples = model.get_samples(100)
+
+        with self.assertRaises(ValueError):
+            engine.simulate(
+                demographic_model=model,
+                contig=contig,
+                samples=samples,
+                dry_run=True,
+            )
+
 
 class TestSlimAvailable(unittest.TestCase):
     """
@@ -536,14 +560,18 @@ def test_engine_available(self):
         self.assertEqual(IS_WINDOWS, "slim" not in all_engines)
 
 
+def get_test_contig(spp="HomSap", chrom="chr22", length_multiplier=0.001):
+    species = stdpopsim.get_species(spp)
+    contig = species.get_contig(chrom, length_multiplier=length_multiplier)
+    contig.fully_neutral(convert_to_substitution=False)
+    return contig
+
+
 class PiecewiseConstantSizeMixin(object):
     """
     Mixin that sets up a simple demographic model used by multiple unit tests.
     """
 
-    species = stdpopsim.get_species("HomSap")
-    contig = species.get_contig("chr22", length_multiplier=0.001)  # ~50 kb
-
     N0 = 1000  # size in the present
     N1 = 500  # ancestral size
     T = 500  # generations since size change occurred
@@ -551,8 +579,8 @@ class PiecewiseConstantSizeMixin(object):
     model = stdpopsim.PiecewiseConstantSize(N0, (T, N1))
     model.generation_time = 1
     samples = model.get_samples(100)
-    mutation_types = [stdpopsim.ext.MutationType(convert_to_substitution=False)]
-    mut_id = len(mutation_types)
+    contig = get_test_contig()
+    mut_id = 0
 
     def allele_frequency(self, ts):
         """
@@ -572,123 +600,189 @@ def allele_frequency(self, ts):
         return af
 
 
+@unittest.skipIf(IS_WINDOWS, "SLiM not available on windows")
+class TestGenomicElementTypes(unittest.TestCase, PiecewiseConstantSizeMixin):
+    def test_default_fully_neutral(self):
+        # if no genomic element type is defined, then by default go fully neutral
+        contig = get_test_contig()
+        contig.genomic_element_types = []
+        contig.mutation_types = []
+        engine = stdpopsim.get_engine("slim")
+        out, _ = capture_output(
+            engine.simulate,
+            demographic_model=self.model,
+            contig=contig,
+            samples=self.samples,
+            slim_script=True,
+        )
+        self.assertEqual(out.count("initializeGenomicElementType"), 1)
+        self.assertEqual(out.count("initializeGenomicElement("), 1)
+
+    def test_single_genomic_element_type_in_script(self):
+        contig = get_test_contig()
+        engine = stdpopsim.get_engine("slim")
+        out, _ = capture_output(
+            engine.simulate,
+            demographic_model=self.model,
+            contig=contig,
+            samples=self.samples,
+            slim_script=True,
+        )
+        self.assertEqual(out.count("initializeGenomicElementType"), 1)
+        self.assertEqual(out.count("initializeGenomicElement("), 1)
+
+    def test_multiple_genomic_element_types_in_script(self):
+        contig = get_test_contig()
+        engine = stdpopsim.get_engine("slim")
+        contig.clear_genomic_mutation_types()
+        contig.add_genomic_element_type(
+            intervals=np.array([[0, 100]]),
+            mutation_types=[stdpopsim.ext.MutationType()],
+            proportions=[0],
+        )
+        contig.add_genomic_element_type(
+            intervals=np.array([[100, 200]]),
+            mutation_types=[stdpopsim.ext.MutationType()],
+            proportions=[0],
+        )
+        contig.add_genomic_element_type(
+            intervals=np.array([[200, 300]]),
+            mutation_types=[stdpopsim.ext.MutationType()],
+            proportions=[0],
+        )
+        out, _ = capture_output(
+            engine.simulate,
+            demographic_model=self.model,
+            contig=contig,
+            samples=self.samples,
+            slim_script=True,
+        )
+        self.assertEqual(out.count("initializeGenomicElementType"), 3)
+        self.assertEqual(out.count("initializeGenomicElement("), 3)
+
+
 @unittest.skipIf(IS_WINDOWS, "SLiM not available on windows")
 class TestMutationTypes(unittest.TestCase, PiecewiseConstantSizeMixin):
     def test_single_mutation_type_in_script(self):
+        contig = get_test_contig()
         engine = stdpopsim.get_engine("slim")
         out, _ = capture_output(
             engine.simulate,
             demographic_model=self.model,
-            contig=self.contig,
+            contig=contig,
             samples=self.samples,
             slim_script=True,
         )
         self.assertEqual(out.count("initializeMutationType"), 1)
 
-        mutation_types = [stdpopsim.ext.MutationType(weight=1)]
         out, _ = capture_output(
             engine.simulate,
             demographic_model=self.model,
-            contig=self.contig,
+            contig=contig,
             samples=self.samples,
-            mutation_types=mutation_types,
             slim_script=True,
         )
         self.assertEqual(out.count("initializeMutationType"), 1)
 
     def test_multiple_mutation_types_in_script(self):
+        contig = get_test_contig()
         engine = stdpopsim.get_engine("slim")
-        mutation_types = [
-            stdpopsim.ext.MutationType(weight=1),
-            stdpopsim.ext.MutationType(weight=2),
+        contig.mutation_types = [
+            stdpopsim.ext.MutationType(),
+            stdpopsim.ext.MutationType(),
         ]
         out, _ = capture_output(
             engine.simulate,
             demographic_model=self.model,
-            contig=self.contig,
+            contig=contig,
             samples=self.samples,
-            mutation_types=mutation_types,
             slim_script=True,
         )
         self.assertEqual(out.count("initializeMutationType"), 2)
 
-        mutation_types = [stdpopsim.ext.MutationType(weight=i) for i in range(10)]
         positive = stdpopsim.ext.MutationType(convert_to_substitution=False)
-        mutation_types.append(positive)
+        contig.mutation_types = [stdpopsim.ext.MutationType() for i in range(10)] + [
+            positive
+        ]
+        contig.genomic_element_types[0].proportions = [1 / 11 for i in range(11)]
+        contig.genomic_element_types[0].mutation_type_ids = [i for i in range(11)]
         out, _ = capture_output(
             engine.simulate,
             demographic_model=self.model,
-            contig=self.contig,
+            contig=contig,
             samples=self.samples,
-            mutation_types=mutation_types,
             slim_script=True,
         )
         self.assertEqual(out.count("initializeMutationType"), 11)
 
     def test_unweighted_mutations_are_not_simulated_by_slim(self):
-        mutation_types = [
+        contig = get_test_contig()
+        contig.mutation_types = [
             stdpopsim.ext.MutationType(convert_to_substitution=True),
             stdpopsim.ext.MutationType(convert_to_substitution=False),
         ]
+        contig.genomic_element_types[0].mutation_type_ids = [0, 1]
+        contig.genomic_element_types[0].proportions = [0, 0]
         ts = slim_simulate_no_recap(
             demographic_model=self.model,
-            contig=self.contig,
+            contig=contig,
             samples=self.samples,
-            mutation_types=mutation_types,
             slim_scaling_factor=10,
             slim_burn_in=0.1,
         )
         self.assertEqual(ts.num_sites, 0)
 
-        mutation_types = [
+        contig.mutation_types = [
             stdpopsim.ext.MutationType(),
             stdpopsim.ext.MutationType(
-                weight=0, distribution_type="g", distribution_args=[-0.01, 0.2]
+                distribution_type="g", distribution_args=[-0.01, 0.2]
             ),
         ]
+        contig.genomic_element_types[0].proportions = [0, 0]
         ts = slim_simulate_no_recap(
             demographic_model=self.model,
-            contig=self.contig,
+            contig=contig,
             samples=self.samples,
-            mutation_types=mutation_types,
             slim_scaling_factor=10,
             slim_burn_in=0.1,
         )
         self.assertEqual(ts.num_sites, 0)
 
     def test_weighted_mutations_are_simulated_by_slim(self):
-        mutation_types = [
-            stdpopsim.ext.MutationType(weight=1, convert_to_substitution=True)
+        contig = get_test_contig()
+        contig.mutation_types = [
+            stdpopsim.ext.MutationType(convert_to_substitution=True)
         ]
+        contig.genomic_element_types[0].proportions = [1]
         ts = slim_simulate_no_recap(
             demographic_model=self.model,
-            contig=self.contig,
+            contig=contig,
             samples=self.samples,
-            mutation_types=mutation_types,
             slim_scaling_factor=10,
             slim_burn_in=0.1,
         )
         self.assertTrue(ts.num_sites > 0)
 
-        mutation_types = [
-            stdpopsim.ext.MutationType(weight=1, convert_to_substitution=False)
+        contig.mutation_types = [
+            stdpopsim.ext.MutationType(convert_to_substitution=False)
         ]
+        contig.genomic_element_types[0].proportions = [1]
         ts = slim_simulate_no_recap(
             demographic_model=self.model,
-            contig=self.contig,
+            contig=contig,
             samples=self.samples,
-            mutation_types=mutation_types,
             slim_scaling_factor=10,
             slim_burn_in=0.1,
         )
         self.assertTrue(ts.num_sites > 0)
 
-        mutation_types = [stdpopsim.ext.MutationType(weight=i) for i in range(10)]
+        contig.mutation_types = [stdpopsim.ext.MutationType() for i in range(10)]
+        contig.genomic_element_types[0].proportions = [1 / 10 for i in range(10)]
+        contig.genomic_element_types[0].mutation_type_ids = list(range(10))
         ts = slim_simulate_no_recap(
             demographic_model=self.model,
-            contig=self.contig,
+            contig=contig,
             samples=self.samples,
-            mutation_types=mutation_types,
             slim_scaling_factor=10,
             slim_burn_in=0.1,
         )
@@ -703,15 +797,6 @@ def test_bad_dominance_coeff(self):
             with self.assertRaises(ValueError):
                 stdpopsim.ext.MutationType(dominance_coeff=dominance_coeff)
 
-    def test_weight(self):
-        for weight in (0, 0.1, 50):
-            stdpopsim.ext.MutationType(weight=weight)
-
-    def test_bad_weight(self):
-        for weight in (-1,):
-            with self.assertRaises(ValueError):
-                stdpopsim.ext.MutationType(weight=weight)
-
     def test_bad_distribution_type(self):
         for distribution_type in (1, {}, None, "~", "!", "F"):
             with self.assertRaises(ValueError):
@@ -743,21 +828,6 @@ def test_bad_distribution_args_g(self):
                     distribution_type="g", distribution_args=distribution_args
                 )
 
-    def test_slim_mutation_frac(self):
-        for weights in ([0.1], [50], [0, 5], [0] * 100 + [0.1]):
-            mut_types = [
-                stdpopsim.ext.MutationType(weight=weight) for weight in weights
-            ]
-            self.assertEqual(1, stdpopsim.ext.slim_mutation_frac(mut_types))
-
-        for weights in ([], [0], [0] * 100):
-            mut_types = [
-                stdpopsim.ext.MutationType(weight=weight) for weight in weights
-            ]
-            self.assertEqual(0, stdpopsim.ext.slim_mutation_frac(mut_types))
-
-        self.assertEqual(0, stdpopsim.ext.slim_mutation_frac(None))
-
 
 @unittest.skipIf(IS_WINDOWS, "SLiM not available on windows")
 class TestDrawMutation(unittest.TestCase, PiecewiseConstantSizeMixin):
@@ -776,7 +846,6 @@ def test_draw_mutation_save(self):
             demographic_model=self.model,
             contig=self.contig,
             samples=self.samples,
-            mutation_types=self.mutation_types,
             extended_events=extended_events,
             dry_run=True,
         )
@@ -795,7 +864,6 @@ def test_draw_mutation_no_save(self):
             demographic_model=self.model,
             contig=self.contig,
             samples=self.samples,
-            mutation_types=self.mutation_types,
             extended_events=extended_events,
             dry_run=True,
         )
@@ -816,7 +884,6 @@ def test_invalid_mutation_type_id(self):
                     demographic_model=self.model,
                     contig=self.contig,
                     samples=self.samples,
-                    mutation_types=self.mutation_types,
                     extended_events=extended_events,
                     dry_run=True,
                 )
@@ -830,11 +897,13 @@ def test_no_mutation_types_defined(self):
                 coordinate=100,
             ),
         ]
+        contig = get_test_contig()
+        contig.mutation_types = []
         engine = stdpopsim.get_engine("slim")
         with self.assertRaises(ValueError):
             engine.simulate(
                 demographic_model=self.model,
-                contig=self.contig,
+                contig=contig,
                 samples=self.samples,
                 extended_events=extended_events,
                 dry_run=True,
@@ -891,7 +960,6 @@ def test_save_point_creation(self):
             demographic_model=self.model,
             contig=self.contig,
             samples=self.samples,
-            mutation_types=self.mutation_types,
             extended_events=extended_events,
             dry_run=True,
         )
@@ -918,7 +986,6 @@ def test_drawn_mutation_not_lost(self):
             demographic_model=self.model,
             contig=self.contig,
             samples=self.samples,
-            mutation_types=self.mutation_types,
             extended_events=extended_events,
             slim_scaling_factor=10,
             slim_burn_in=0.1,
@@ -947,7 +1014,6 @@ def test_drawn_mutation_is_lost(self):
             demographic_model=self.model,
             contig=self.contig,
             samples=self.samples,
-            mutation_types=self.mutation_types,
             extended_events=extended_events,
             slim_scaling_factor=10,
             slim_burn_in=0.1,
@@ -978,7 +1044,6 @@ def test_drawn_mutation_meets_AF_threshold(self):
                 demographic_model=self.model,
                 contig=self.contig,
                 samples=self.samples,
-                mutation_types=self.mutation_types,
                 extended_events=extended_events,
                 slim_scaling_factor=10,
                 slim_burn_in=0.1,
@@ -1062,7 +1127,6 @@ def test_bad_GenerationAfter_times(self):
                     demographic_model=self.model,
                     contig=self.contig,
                     samples=self.samples,
-                    mutation_types=self.mutation_types,
                     extended_events=extended_events,
                     dry_run=True,
                 )
@@ -1100,7 +1164,6 @@ def test_conditioning_without_save(self):
                 demographic_model=self.model,
                 contig=self.contig,
                 samples=self.samples,
-                mutation_types=self.mutation_types,
                 extended_events=extended_events,
                 slim_scaling_factor=10,
                 slim_burn_in=0.1,
@@ -1124,7 +1187,6 @@ def test_no_drawn_mutation(self):
                 demographic_model=self.model,
                 contig=self.contig,
                 samples=self.samples,
-                mutation_types=self.mutation_types,
                 extended_events=extended_events,
                 dry_run=True,
             )
@@ -1181,7 +1243,6 @@ def test_positive_mutation_meets_AF_threshold(self):
                 demographic_model=self.model,
                 contig=self.contig,
                 samples=self.samples,
-                mutation_types=self.mutation_types,
                 extended_events=extended_events,
                 slim_scaling_factor=10,
                 slim_burn_in=0.1,
@@ -1207,7 +1268,6 @@ def test_no_drawn_mutation(self):
                 demographic_model=self.model,
                 contig=self.contig,
                 samples=self.samples,
-                mutation_types=self.mutation_types,
                 extended_events=extended_events,
                 dry_run=True,
             )
@@ -1259,7 +1319,6 @@ def test_bad_GenerationAfter_times(self):
                     demographic_model=self.model,
                     contig=self.contig,
                     samples=self.samples,
-                    mutation_types=self.mutation_types,
                     extended_events=extended_events,
                     dry_run=True,
                 )
@@ -1283,3 +1342,138 @@ def test_bad_extended_events(self):
                     extended_events=[bad_ee],
                     dry_run=True,
                 )
+
+
+class TestMiscFunctions(unittest.TestCase):
+    def get_test_contig(self):
+        species = stdpopsim.get_species("HomSap")
+        contig = species.get_contig("chr22", length_multiplier=0.001)
+        return contig
+
+    def test_slim_fractions(self):
+        contig = self.get_test_contig()
+
+        self.assertTrue(stdpopsim.slim_engine.get_slim_fractions(contig).size == 0)
+        contig.fully_neutral()
+        self.assertTrue(
+            np.isclose(stdpopsim.slim_engine.get_slim_fractions(contig), np.array([0]))
+        )
+        contig.fully_neutral(slim_mutations=True)
+        self.assertTrue(
+            np.isclose(stdpopsim.slim_engine.get_slim_fractions(contig), np.array([1]))
+        )
+
+        proportions = ([0, 1.0], [0.5, 0.5], [0, 0], [0.2, 0])
+        for props in proportions:
+            slim_frac = np.array(sum(props))
+            contig.clear_genomic_mutation_types()
+            contig.add_genomic_element_type(
+                intervals=np.array([[0, 1]]),
+                mutation_types=[
+                    stdpopsim.ext.MutationType(),
+                    stdpopsim.ext.MutationType(),
+                ],
+                proportions=props,
+            )
+            self.assertTrue(
+                (slim_frac == stdpopsim.slim_engine.get_slim_fractions(contig)).all()
+            )
+
+    def test_msp_mutation_rate_map(self):
+        contig = self.get_test_contig()
+        contig.fully_neutral()
+        rmap, _ = stdpopsim.slim_engine.get_msp_and_slim_mutation_rate_maps(contig)
+        self.assertTrue(
+            np.allclose(
+                rmap.position,
+                np.array([0, int(contig.recombination_map.sequence_length)]),
+            )
+        )
+        self.assertTrue(np.allclose(rmap.rate, np.array([contig.mutation_rate])))
+
+    def test_slim_mutation_rate_map(self):
+        contig = self.get_test_contig()
+        contig.fully_neutral()
+        _, (breaks, rates) = stdpopsim.slim_engine.get_msp_and_slim_mutation_rate_maps(
+            contig
+        )
+        self.assertTrue(breaks == [int(contig.recombination_map.sequence_length) - 1])
+        self.assertTrue(rates == [0.0])
+
+    def test_complex_mutation_rate_maps(self):
+        contig = self.get_test_contig()
+        for prop1, prop2 in ((0.3, 0.1), (0.7, 0.2), (1.0, 0.0)):
+            contig.clear_genomic_mutation_types()
+            contig.add_genomic_element_type(
+                intervals=np.array([[10, 100]]),
+                mutation_types=[stdpopsim.ext.MutationType()],
+                proportions=[prop1],
+            )
+
+            rmap, (
+                obs_slim_breaks,
+                obs_slim_rates,
+            ) = stdpopsim.slim_engine.get_msp_and_slim_mutation_rate_maps(contig)
+            obs_msp_breaks, obs_msp_rates = rmap.position, rmap.rate
+
+            exp_msp_breaks = [0, 10, 100, int(contig.recombination_map.sequence_length)]
+            exp_msp_rates = [
+                contig.mutation_rate,
+                contig.mutation_rate * (1 - prop1),
+                contig.mutation_rate,
+            ]
+            self.assertTrue(np.allclose(obs_msp_breaks, exp_msp_breaks))
+            self.assertTrue(
+                (np.isclose(np.array(obs_msp_rates), np.array(exp_msp_rates))).all()
+            )
+
+            exp_slim_breaks = [9, 99, int(contig.recombination_map.sequence_length) - 1]
+            exp_slim_rates = [0, contig.mutation_rate * prop1, 0]
+            self.assertEqual(obs_slim_breaks, exp_slim_breaks)
+            self.assertTrue(
+                (np.isclose(np.array(obs_slim_rates), np.array(exp_slim_rates))).all()
+            )
+
+            contig.clear_genomic_mutation_types()
+            contig.add_genomic_element_type(
+                intervals=np.array([[0, 50]]),
+                mutation_types=[stdpopsim.ext.MutationType()],
+                proportions=[prop1],
+            )
+            contig.add_genomic_element_type(
+                intervals=np.array([[50, 100]]),
+                mutation_types=[stdpopsim.ext.MutationType()],
+                proportions=[prop2],
+            )
+
+            rmap, (
+                obs_slim_breaks,
+                obs_slim_rates,
+            ) = stdpopsim.slim_engine.get_msp_and_slim_mutation_rate_maps(contig)
+            obs_msp_breaks, obs_msp_rates = rmap.position, rmap.rate
+
+            exp_msp_breaks = [0, 50, 100, int(contig.recombination_map.sequence_length)]
+            exp_msp_rates = [
+                contig.mutation_rate * (1 - prop1),
+                contig.mutation_rate * (1 - prop2),
+                contig.mutation_rate,
+            ]
+            self.assertTrue(np.allclose(obs_msp_breaks, exp_msp_breaks))
+            self.assertTrue(
+                (np.isclose(np.array(obs_msp_rates), np.array(exp_msp_rates))).all()
+            )
+
+            exp_slim_breaks = [
+                49,
+                99,
+                int(contig.recombination_map.sequence_length) - 1,
+            ]
+            exp_slim_rates = [
+                contig.mutation_rate * prop1,
+                contig.mutation_rate * prop2,
+                0,
+            ]
+            self.assertEqual(obs_slim_breaks, exp_slim_breaks)
+            self.assertTrue(
+                (np.isclose(np.array(obs_slim_rates), np.array(exp_slim_rates))).all()
+            )
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 4c8b7c701..d2ed2580e 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -6,6 +6,7 @@
 import pathlib
 import tarfile
 import tempfile
+import numpy as np
 
 from stdpopsim import utils
 from stdpopsim import Chromosome, Genome
@@ -434,3 +435,62 @@ def test_add_only_unique_synonyms(self):
         utils.append_common_synonyms(genome)
         self.assertTrue("chr1" in genome.chromosomes[0].synonyms)
         self.assertEqual(len(genome.chromosomes[0].synonyms), 1)
+
+
+class TestIntervalUtilities(unittest.TestCase):
+    def test_intervals_array_shape(self):
+        # shape (n, x) where x >= 2
+        bad_arrays = (np.array([10, 20]), np.array([[[10, 20]]]))
+        for bad_array in bad_arrays:
+            with self.assertRaises(ValueError):
+                utils.check_intervals_array_shape(intervals=bad_array)
+            with self.assertRaises(ValueError):
+                utils.build_intervals_array(intervals=bad_array)
+        good_arrays = (
+            np.array([[10, 20], [30, 40]]),
+            np.array([[10, 20, 1], [30, 40, 2]]),
+        )
+        for good_array in good_arrays:
+            utils.check_intervals_array_shape(intervals=good_array)
+            built_array = utils.build_intervals_array(intervals=good_array)
+            self.assertTrue((built_array == good_array).all())
+
+    def test_invalid_intervals(self):
+        # right <= left, overlapping, uncastable float->int
+        invalid_arrays = (
+            np.array([[20, 10]]),
+            np.array([[10, 20], [15, 20]]),
+            np.array([[10, 20], [13, 41]]),
+        )
+        for invalid_array in invalid_arrays:
+            with self.assertRaises(ValueError):
+                utils.build_intervals_array(intervals=invalid_array)
+            with self.assertRaises(ValueError):
+                utils.check_intervals_validity(intervals=invalid_array)
+        # interval outside [start, end)
+        with self.assertRaises(ValueError):
+            utils.build_intervals_array(intervals=[[10, 50]], start=20)
+        with self.assertRaises(ValueError):
+            utils.check_intervals_validity(intervals=np.array([[10, 50]]), end=30)
+
+    def test_intervals_casting(self):
+        castable_intervals = (
+            [[10, 20], [20, 40]],
+            [[1, 2, 2], [2, 3, 2]],
+            np.array([[10, 20]]),
+        )
+        for intervals in castable_intervals:
+            casted = utils.build_intervals_array(intervals)
+            self.assertTrue(isinstance(casted, np.ndarray) and casted.dtype == np.int64)
+            self.assertTrue(casted.shape[0] == len(intervals))
+            self.assertTrue(casted.shape[1] >= 2)
+
+    def test_interval_sorting(self):
+        unsorted_intervals = (
+            np.array([[20, 30], [10, 20]]),
+            np.array([[10, 20], [20, 30], [2, 3]]),
+        )
+        for intervals in unsorted_intervals:
+            casted = utils.build_intervals_array(intervals)
+            self.assertFalse(np.all(np.diff(intervals[:, 0]) >= 0))
+            self.assertTrue(np.all(np.diff(casted[:, 0]) >= 0))