diff --git a/atomica/demos.py b/atomica/demos.py
index b56ed912..c9e30de0 100644
--- a/atomica/demos.py
+++ b/atomica/demos.py
@@ -11,10 +11,30 @@
 from .scenarios import BudgetScenario
 from .utils import TimeSeries
 
-__all__ = ["demo", "make_demo_scenarios"]
-
-
-def demo(which: str = None, do_run: bool = True, addprogs: bool = True) -> Project:
+__all__ = ["demo", "demo_projects", "make_demo_scenarios"]
+
+
+demo_projects = [
+    "udt",
+    "udt_dyn",
+    "usdt",
+    "cervicalcancer",
+    "sir",
+    "diabetes",
+    "combined",
+    "service",
+    "hypertension",
+    "hypertension_dyn",
+    "hiv",
+    "hiv_dyn",
+    "tb_simple",
+    "tb_simple_dyn",
+    "tb",
+    "dt",
+]
+
+
+def demo(which: str = "sir", do_run: bool = True, addprogs: bool = True) -> Project:
     """
     Return a demo project
 
@@ -25,30 +45,11 @@ def demo(which: str = None, do_run: bool = True, addprogs: bool = True) -> Proje
 
     """
 
-    options = [
-        "udt",
-        "udt_dyn",
-        "usdt",
-        "cervicalcancer",
-        "sir",
-        "diabetes",
-        "combined",
-        # 'service',
-        "hypertension",
-        "hypertension_dyn",
-        "hiv",
-        "hiv_dyn",
-        "tb_simple",
-        "tb_simple_dyn",
-        "environment",
-        "tb",
-    ]
-
-    dtdict = sc.odict.fromkeys(options, 1.0)
+    dtdict = sc.odict.fromkeys(demo_projects, 1.0)
     dtdict["tb"] = 0.5
 
-    if which is None or which not in options:
-        raise Exception("Supported project types are:\n%s" % ("\n".join(options)))
+    if which not in demo_projects:
+        raise Exception("Supported demonstration projects are:\n%s" % ("\n".join(demos)))
 
     framework = LIBRARY_PATH / f"{which}_framework.xlsx"
     databook = LIBRARY_PATH / f"{which}_databook.xlsx"
@@ -61,7 +62,7 @@ def demo(which: str = None, do_run: bool = True, addprogs: bool = True) -> Proje
     if do_run:
         P.run_sim(store_results=True)
 
-    if addprogs:
+    if addprogs and progbook.exists():
         logger.debug("Loading progbook")
         P.load_progbook(progbook)
 
diff --git a/atomica/excel.py b/atomica/excel.py
index d8b63c98..58c92792 100644
--- a/atomica/excel.py
+++ b/atomica/excel.py
@@ -448,7 +448,8 @@ def __init__(self, code_name: str, full_name: str, tvec: np.array, from_pops: li
             raise Exception('Unknown TimeDependentConnections type - must be "transfer" or "interaction"')
 
     def __repr__(self):
-        return '<TDC %s "%s">' % (self.type.title(), self.code_name)
+        # return '<TDC %s "%s">' % (self.type.title(), self.code_name)
+        return sc.prepr(self)
 
     @classmethod
     def from_tables(cls, tables: list, interaction_type):
@@ -862,7 +863,7 @@ def _write_pop_matrix(self, worksheet, start_row, formats, references: dict = No
             from_idx = self.from_pops.index(from_pop)
             to_idx = self.to_pops.index(to_pop)
             if boolean_choice:
-                value = "Y" if value else "N"
+                value = "Y" if (value.vals or value.assumption) else "N"
             content[from_idx, to_idx] = value
 
         # Write the content
diff --git a/atomica/migration.py b/atomica/migration.py
index 4d6c52a9..673dba3e 100644
--- a/atomica/migration.py
+++ b/atomica/migration.py
@@ -868,4 +868,5 @@ def _progset_nontargetable_flag(progset):
     for d in progset.comps.values():
         if 'non_targetable' not in d:
             d['non_targetable'] = False
-    return progset
\ No newline at end of file
+    return progset
+
diff --git a/atomica/model.py b/atomica/model.py
index c89b6797..2a10425f 100644
--- a/atomica/model.py
+++ b/atomica/model.py
@@ -22,6 +22,8 @@
 from .parameters import Parameter as ParsetParameter
 from .parameters import ParameterSet as ParameterSet
 import math
+
+from .utils import stochastic_rounding
 import pandas as pd
 
 model_settings = dict()
@@ -229,11 +231,15 @@ def __setitem__(self, key, value) -> None:
 class Compartment(Variable):
     """A class to wrap up data for one compartment within a cascade network."""
 
-    def __init__(self, pop, name):
+    def __init__(self, pop, name, multinomial: bool = False, rng_sampler=None):
         Variable.__init__(self, pop=pop, id=(pop.name, name))
         self.units = "Number of people"
         self.outlinks = []
         self.inlinks = []
+        self.multinomial = multinomial
+        self.rng_sampler = rng_sampler
+        if self.multinomial and rng_sampler is None:
+            self.rng_sampler = np.random.default_rng()
         self._cached_outflow = None
 
     def unlink(self):
@@ -353,11 +359,21 @@ def resolve_outflows(self, ti: int) -> None:
         else:
             rescale = 1
 
-        n = rescale * self.vals[ti]
         self._cached_outflow = 0
-        for link in self.outlinks:
-            link.vals[ti] = link._cache * n
-            self._cached_outflow += link.vals[ti]
+
+        if self.multinomial:
+            # for every person in the compartment, we now have probabilities for each possible outflow or staying as the remainder
+            stays = max(0.0, 1.0 - outflow)
+            rescaled_frac_flows = [rescale * link._cache for link in self.outlinks] + [stays]
+            number_outflows = self.rng_sampler.multinomial(self[ti], rescaled_frac_flows, size=1)[0]
+            for ln, link in enumerate(self.outlinks):
+                link.vals[ti] = number_outflows[ln]
+                self._cached_outflow += link.vals[ti]
+        else:
+            n = rescale * self.vals[ti]
+            for link in self.outlinks:
+                link.vals[ti] = link._cache * n
+                self._cached_outflow += link.vals[ti]
 
     def update(self, ti: int) -> None:
         """
@@ -399,7 +415,7 @@ def connect(self, dest, par) -> None:
 
 
 class JunctionCompartment(Compartment):
-    def __init__(self, pop, name: str, duration_group: str = None):
+    def __init__(self, pop, name: str, duration_group: str = None, multinomial: bool = False, rng_sampler=None):
         """
 
         A TimedCompartment has a duration group by virtue of having a `.parameter` attribute and a flush link.
@@ -418,7 +434,7 @@ def __init__(self, pop, name: str, duration_group: str = None):
         :param duration_group: Optionally specify a duration group
 
         """
-        super().__init__(pop, name)
+        super().__init__(pop, name, multinomial=multinomial, rng_sampler=rng_sampler)
         self.duration_group = duration_group  #: Store the name of the duration group, if the junction belongs to one
 
     def connect(self, dest, par) -> None:
@@ -517,12 +533,21 @@ def balance(self, ti: int) -> None:
         outflow_fractions = [link.parameter.vals[ti] for link in self.outlinks]
         total_outflow = sum(outflow_fractions)
 
-        # Finally, assign the inflow to the outflow proportionately accounting for the total outflow downscaling
+        outflow_fractions /= total_outflow  # proportionately accounting for the total outflow downscaling
+
+        # assign outflow stochastically to compartments based on the probabilities
+        if self.multinomial and net_inflow > 0.0:
+            if self.duration_group:
+                raise Exception("Need to implement code here to make sure outflows from EACH part of the duration group are integers")
+            else:
+                outflow_fractions = self.rng_sampler.multinomial(net_inflow, outflow_fractions, size=1)[0] / net_inflow
+
+        # Finally, assign the inflow to the outflow proportionately
         for frac, link in zip(outflow_fractions, self.outlinks):
             if self.duration_group:
-                link._vals[:, ti] = net_inflow * frac / total_outflow
+                link._vals[:, ti] = net_inflow * frac
             else:
-                link.vals[ti] = net_inflow * frac / total_outflow
+                link.vals[ti] = net_inflow * frac
 
     def initial_flush(self) -> None:
         """
@@ -543,6 +568,10 @@ def initial_flush(self) -> None:
             outflow_fractions = np.array([link.parameter.vals[0] for link in self.outlinks])
             outflow_fractions /= np.sum(outflow_fractions)
 
+            # assign outflow stochastically to compartments based on the probabilities
+            if self.multinomial:
+                outflow_fractions = self.rng_sampler.multinomial(self.vals[0], outflow_fractions, size=1)[0] / self.vals[0]
+
             # Assign the inflow directly to the outflow compartments
             # This is done using the [] indexing on the downstream compartment so it is agnostic
             # to the downstream compartment type
@@ -602,22 +631,37 @@ def balance(self, ti: int) -> None:
             net_inflow = net_inflow + link._vals[:, ti]  # Get the flow from TimedLink._vals. nb. using += doesn't work with some array size combinations
 
         outflow_fractions = np.zeros(len(self.outlinks))
+        residual_link_index = None
         for i, link in enumerate(self.outlinks):
             if link.parameter is not None:
                 outflow_fractions[i] = link.parameter.vals[ti]
             else:
+                residual_link_index = i
                 outflow_fractions[i] = 0
 
         total_outflow = outflow_fractions.sum()
         if total_outflow > 1:
             outflow_fractions /= total_outflow
 
+        # assign outflow stochastically to compartments based on the probabilities
+        if self.multinomial and net_inflow > 0.0:
+            if np.abs(np.round(net_inflow) - net_inflow).sum() > 1e-3:
+                raise Exception(f"Net inflow should be an integer for a stochastic run, was: {net_inflow}")
+
+            # Give residual probability to the residual link if there is one, otherwise normalize probabilities
+            if residual_link_index is not None:
+                outflow_fractions[residual_link_index] = max(0, 1 - outflow_fractions.sum())
+            else:  # Probabilities have to sum to one as we have to flush the junction, and multinomials probabilities should sum to 1
+                outflow_fractions = outflow_fractions / outflow_fractions.sum()
+
+            outflow_fractions = self.rng_sampler.multinomial(np.round(net_inflow).astype(int), outflow_fractions, size=1)[0] / (net_inflow)
+
         # Calculate outflows
         outflow = net_inflow.reshape(-1, 1) * outflow_fractions.reshape(1, -1)
 
         # Finally, assign the inflow to the outflow proportionately accounting for the total outflow downscaling
         for frac, link in zip(outflow_fractions, self.outlinks):
-            if link.parameter is None and total_outflow < 1:
+            if link.parameter is None and total_outflow < 1 and not self.multinomial:  # If we are stochastic then we have already taken the residual link into account
                 flow = net_inflow - np.sum(outflow, axis=1)  # Sum after multiplying by outflow fractions to reduce numerical precision errors and enforce conserved quantities more accurately
             else:
                 flow = net_inflow * frac
@@ -647,6 +691,10 @@ def initial_flush(self) -> None:
                 outflow_fractions /= total_outflow
                 has_residual = False
 
+            # assign outflow stochastically to compartments based on the probabilities
+            if self.multinomial and total_outflow > 0.0:
+                outflow_fractions = self.rng_sampler.multinomial(self.vals[0], outflow_fractions, size=1)[0] / (self.vals[0] * sum(outflow_fractions))
+
             # Assign the inflow directly to the outflow compartments
             # This is done using the [] indexing on the downstream compartment so it is agnostic
             # to the downstream compartment type
@@ -677,24 +725,28 @@ class SourceCompartment(Compartment):
 
     """
 
-    def __init__(self, pop, name):
-        super().__init__(pop, name)
+    def __init__(self, pop, name, multinomial: bool = False, rng_sampler=None):
+        super().__init__(pop, name, multinomial=multinomial, rng_sampler=rng_sampler)
 
     def preallocate(self, tvec: np.array, dt: float) -> None:
         super().preallocate(tvec, dt)
-        self.vals.fill(0.0)  #: Source compartments have an unlimited number of people in them. TODO: If this complicates validation, set to 0.0
+        self.vals.fill(0.0)  #: Source compartments have an unlimited number of people in them. Set to 0.0 to simplify validation.
 
     def resolve_outflows(self, ti: int) -> None:
+        # Unlike other Compartments, link._cache will still be in Number form and not converted to a proportion
         for link in self.outlinks:
-            link.vals[ti] = link._cache
+            if self.multinomial:
+                link.vals[ti] = stochastic_rounding(link._cache, rng_sampler=self.rng_sampler)
+            else:
+                link.vals[ti] = link._cache
 
     def update(self, ti: int) -> None:
         pass
 
 
 class SinkCompartment(Compartment):
-    def __init__(self, pop, name):
-        super().__init__(pop, name)
+    def __init__(self, pop, name, multinomial: bool = False, rng_sampler=None):
+        super().__init__(pop, name, multinomial=multinomial, rng_sampler=rng_sampler)
 
     def preallocate(self, tvec: np.array, dt: float) -> None:
         """
@@ -750,7 +802,7 @@ def update(self, ti: int) -> None:
 
 
 class TimedCompartment(Compartment):
-    def __init__(self, pop, name: str, parameter: ParsetParameter):
+    def __init__(self, pop, name: str, parameter: ParsetParameter, multinomial: bool = False, rng_sampler=None):
         """
         Instantiate the TimedCompartment
 
@@ -763,7 +815,7 @@ def __init__(self, pop, name: str, parameter: ParsetParameter):
 
         """
 
-        Compartment.__init__(self, pop=pop, name=name)
+        Compartment.__init__(self, pop=pop, name=name, multinomial=multinomial, rng_sampler=rng_sampler)
         self._vals = None  #: Primary storage, a matrix of size (duration x timesteps). The first row is the one that people leave from, the last row is the one that people arrive in
         self.parameter = parameter  #: The parameter to read the duration from - this needs to be done after parameters are precomputed though, in case the duration is coming from a (constant) function
         self.flush_link = None  #: Reference to the timed outflow link that flushes the compartment. Note that this needs to be set externally because Links are instantiated after Compartments
@@ -867,6 +919,9 @@ def resolve_outflows(self, ti: int) -> None:
 
         """
 
+        # TODO adjust for stochastic variables
+        assert self.multinomial == False, "resolve_outflows needs to be updated to work with stochastic transitions"
+
         # First, work out the scale factors as usual
         self.flush_link._cache = 0.0  # At this stage, no outflow goes via the flush link
 
@@ -1646,7 +1701,7 @@ class Population:
     Each model population must contain a set of compartments with equivalent names.
     """
 
-    def __init__(self, framework, name: str, label: str, progset: ProgramSet, pop_type: str, build_cache: dict | None = None):
+    def __init__(self, framework, name: str, label: str, progset: ProgramSet, pop_type: str, multinomial: bool = False, rng_sampler=None, build_cache: dict | None = None):
         """
         Construct a Population
 
@@ -1660,11 +1715,16 @@ def __init__(self, framework, name: str, label: str, progset: ProgramSet, pop_ty
         :param label: The full name of the population
         :param progset: A ProgramSet instance
 
+        :param multinomial: whether to use stochasticity in determining population movement
+        :param rng_sampler: optional Generator for random numbers to give consistent results between model runs
+
         """
 
         self.name = name  #: The code name of the population
         self.label = label  #: The full name/label of the population
         self.type = pop_type  #: The population's type
+        self.multinomial = multinomial  #: Whether the population should be handled stochastically as discrete integer people instead of fractions
+        self.rng_sampler = rng_sampler
 
         self.comps = list()  #: List of Compartment objects
         self.characs = list()  #: List of Characteristic objects
@@ -1882,17 +1942,17 @@ def build(self, framework, progset, build_cache: dict | None = None):
             population_type, duration_group, is_junction, is_source, is_sink = comp_rows[comp_name]
             if population_type == self.type:
                 if comp_name in residual_junctions:
-                    self.comps.append(ResidualJunctionCompartment(pop=self, name=comp_name, duration_group=duration_group))
+                    self.comps.append(ResidualJunctionCompartment(pop=self, name=comp_name, multinomial=self.multinomial, rng_sampler=self.rng_sampler, duration_group=duration_group))
                 elif is_junction == "y":
-                    self.comps.append(JunctionCompartment(pop=self, name=comp_name, duration_group=duration_group))
+                    self.comps.append(JunctionCompartment(pop=self, name=comp_name, multinomial=self.multinomial, rng_sampler=self.rng_sampler, duration_group=duration_group))
                 elif duration_group:
-                    self.comps.append(TimedCompartment(pop=self, name=comp_name, parameter=self.par_lookup[duration_group]))
+                    self.comps.append(TimedCompartment(pop=self, name=comp_name, multinomial=self.multinomial, rng_sampler=self.rng_sampler, parameter=self.par_lookup[duration_group]))
                 elif is_source == "y":
-                    self.comps.append(SourceCompartment(pop=self, name=comp_name))
+                    self.comps.append(SourceCompartment(pop=self, name=comp_name, multinomial=self.multinomial, rng_sampler=self.rng_sampler))
                 elif is_sink == "y":
-                    self.comps.append(SinkCompartment(pop=self, name=comp_name))
+                    self.comps.append(SinkCompartment(pop=self, name=comp_name, multinomial=self.multinomial, rng_sampler=self.rng_sampler))
                 else:
-                    self.comps.append(Compartment(pop=self, name=comp_name))
+                    self.comps.append(Compartment(pop=self, name=comp_name, multinomial=self.multinomial, rng_sampler=self.rng_sampler))
 
         self.comp_lookup = {comp.name: comp for comp in self.comps}
 
@@ -2122,13 +2182,16 @@ def report_characteristic(charac, n_indent=0):
 
         # Insert the calculated initial values
         for i, c in enumerate(comps):
-            c[0] = max(0.0, x[i,0])
+            if self.multinomial:
+                c[0] = stochastic_rounding(max(0.0, x[i,0]), rng_sampler=self.rng_sampler)  # integer values (still represented as floats)
+            else:
+                c[0] = max(0.0, x[i,0])
 
 
 class Model:
     """A class to wrap up multiple populations within model and handle cross-population transitions."""
 
-    def __init__(self, settings, framework, parset, progset=None, program_instructions=None):
+    def __init__(self, settings, framework, parset, progset=None, program_instructions=None, rng_sampler=None, acceptance_criteria=[]):
 
         # Note that if a progset is provided and program instructions are not, then programs will not be
         # turned on. However, the progset is still available so that the coverage can still be probed
@@ -2148,6 +2211,13 @@ def __init__(self, settings, framework, parset, progset=None, program_instructio
         self.program_instructions = sc.dcp(program_instructions)  # program instructions
         self.t = settings.tvec  #: Simulation time vector (this is a brand new instance from the `settings.tvec` property method)
         self.dt = settings.sim_dt  #: Simulation time step
+        self.multinomial = settings.multinomial  #: Run with discrete numbers of people per compartment instead of fractions.
+
+        self.acceptance_criteria = sc.dcp(acceptance_criteria)  #: If the model fails to adhere to these criteria, raise a BadInitialization exception
+
+        if rng_sampler is None:
+            rng_sampler = np.random.default_rng()
+        self.rng_sampler = rng_sampler
 
         self._t_index = 0  # Keeps track of array index for current timepoint data within all compartments.
         self._vars_by_pop = None  # Cache to look up lists of variables by name across populations: {var_name:[var instances]}
@@ -2295,7 +2365,7 @@ def build(self, parset):
         # First construct populations
         build_cache = _get_framework_build_cache(self.framework)
         for k, (pop_name, pop_label, pop_type) in enumerate(zip(parset.pop_names, parset.pop_labels, parset.pop_types)):
-            self.pops.append(Population(framework=self.framework, name=pop_name, label=pop_label, progset=self.progset, pop_type=pop_type, build_cache=build_cache))
+            self.pops.append(Population(framework=self.framework, name=pop_name, label=pop_label, progset=self.progset, pop_type=pop_type, multinomial=self.multinomial, rng_sampler=self.rng_sampler, build_cache=build_cache))
             self._pop_ids[pop_name] = k
 
         # Expand interactions into matrix form
@@ -2419,6 +2489,21 @@ def build(self, parset):
                 obj.preallocate(self.t, self.dt)
             pop.initialize_compartments(parset, self.framework, self.t[0])
 
+        # More finally, set up acceptance criteria for those which should be checked at runtime (variables that have ._is_dynamic = True), and those which can only be checked after conclusion
+        if self.acceptance_criteria:
+            for criteria in self.acceptance_criteria:
+                pars = self._vars_by_pop[criteria["parameter"]]
+                for par in pars:
+                    if par.pop.name == criteria["population"]:
+                        if par._is_dynamic:
+                            criteria["evaluate_at"] = self.t[np.where(self.t > criteria["t_range"][1])][0]
+                        else:
+                            criteria["evaluate_at"] = np.inf  # can only evaluate postcompute
+
+            # sort by evaluation time - now we only have to look at the first element every timestep
+            self.acceptance_criteria = sorted(self.acceptance_criteria, key=lambda item: item["evaluate_at"])
+            self.acceptance_index = 0  # which index to evaluate
+
     def _set_exec_order(self) -> None:
         """
         Get the execution order
@@ -2542,6 +2627,8 @@ def process(self) -> None:
             self.update_comps()
             self.update_pars()
             self.update_links()
+            if self.acceptance_criteria:
+                self.update_acceptance()
 
         # Update postcompute parameters - note that it needs to be done in execution order
         for par_name in self._exec_order["all_pars"]:
@@ -2550,6 +2637,9 @@ def process(self) -> None:
                     par.update()
                     par.constrain()
 
+        if self.acceptance_criteria:  # call update_acceptance once more for any variables that could only be assessed postcompute
+            self.update_acceptance(evaluate_all=True)
+
         # Clear characteristic internal storage and switch to dynamic computation to save space
         for pop in self.pops:
             for charac in pop.characs:
@@ -2677,6 +2767,39 @@ def update_comps(self) -> None:
             for comp in pop.comps:
                 comp.update(ti)
 
+    def update_acceptance(self, evaluate_all=False) -> None:
+        """
+        Update any acceptance criteria for the model run, raise BadInitialization error if the model run is failing
+
+        :param evaluate_all: Final evaluation of all acceptance criteria.
+
+        """
+
+        time = self.t[self._t_index]  # actual time
+
+        # We know that these are sorted by evaluation time so we only need to look at the first index
+        while (self.acceptance_index < len(self.acceptance_criteria)) and (evaluate_all or time >= self.acceptance_criteria[self.acceptance_index]["evaluate_at"]):
+            # time ranges for assessment are in order, so only need to check the first one
+            assessable = self.acceptance_criteria[self.acceptance_index]
+            # we have passed the time window, assess if the criteria was met - either return a fail condition or remove this acceptance criteria
+            a_par = assessable["parameter"]
+            a_pop = assessable["population"]
+            a_times = assessable["t_range"]
+            accept_vals = assessable["value"]
+            a_t_inds = np.where(np.logical_and(self.t >= a_times[0], self.t <= a_times[1]))  # self.t is the full tvec for the model
+
+            pars = self._vars_by_pop[a_par]
+            for par in pars:  # TODO should be more efficient way to directly access the par values for the pop rather than looping
+                if par.pop.name == a_pop:
+                    a_t_val = np.mean(par[a_t_inds])  # mean across the time range
+                    if a_t_val <= accept_vals[0] or a_t_val >= accept_vals[1]:
+                        # print(f'{time} Rejecting run as sampled sim value {a_t_val} outside of acceptable bound {accept_vals} for parameter {a_par}, population {a_pop} at time {a_times}')
+                        raise BadInitialization(f"Rejecting run as sampled sim value {a_t_val} outside of acceptable bound {accept_vals} for parameter {a_par}, population {a_pop} at time {a_times}")
+                    # else:
+                    # print(f'{time} ACCEPTING run as sampled sim value {a_t_val} inside of acceptable bound {accept_vals} for parameter {a_par}, population {a_pop} at time {a_times}')
+            assessable["assessed_value"] = a_t_val
+            self.acceptance_index += 1
+
     def flush_junctions(self) -> None:
         """
         Flush initialization values from junctions
@@ -2850,7 +2973,7 @@ def update_pars(self) -> None:
                     par.constrain(ti)
 
 
-def run_model(settings, framework, parset: ParameterSet, progset: ProgramSet = None, program_instructions: ProgramInstructions = None, name: str = None):
+def run_model(settings, framework, parset: ParameterSet, progset: ProgramSet = None, program_instructions: ProgramInstructions = None, name: str = None, rng_sampler=None, acceptance_criteria=[]):
     """
     Build and process model
 
@@ -2874,10 +2997,12 @@ def run_model(settings, framework, parset: ParameterSet, progset: ProgramSet = N
     :param progset: Optionally provide a :class:`ProgramSet` instance to use programs
     :param program_instructions: Optional :class:`ProgramInstructions` instance. If ``progset`` is specified, then instructions must be provided
     :param name: Optionally specify the name to assign to the output result
+    :param rng_sampler: Optionally specify a random number generator that may have been seeded to generate consistent results
+
     :return: A :class:`Result` object containing the processed model
 
     """
 
-    m = Model(settings, framework, parset, progset, program_instructions)
+    m = Model(settings, framework, parset, progset, program_instructions, rng_sampler=rng_sampler, acceptance_criteria=acceptance_criteria)
     m.process()
     return Result(model=m, parset=parset, name=name)
diff --git a/atomica/parameters.py b/atomica/parameters.py
index 4f88b15f..7bf32e89 100644
--- a/atomica/parameters.py
+++ b/atomica/parameters.py
@@ -96,7 +96,7 @@ def interpolate(self, tvec, pop_name: str) -> np.array:
 
         return self.ts[pop_name].interpolate(tvec, method=self._interpolation_method)
 
-    def sample(self, constant: bool) -> None:
+    def sample(self, constant: bool, rng_sampler=None) -> None:
         """
         Perturb parameter based on uncertainties
 
@@ -105,11 +105,11 @@ def sample(self, constant: bool) -> None:
 
         :param constant: If True, time series will be perturbed by a single constant offset. If False,
                          a different perturbation will be applied to each time specific value independently.
-
+        :param rng_sampler: Optional random number generator that may have been seeded to generate consistent results
         """
 
         for k, ts in self.ts.items():
-            self.ts[k] = ts.sample(constant)
+            self.ts[k] = ts.sample(constant, rng_sampler)
 
     def smooth(self, tvec, method="smoothinterp", pop_names=None, **kwargs):
         """
@@ -530,19 +530,20 @@ def get_par(self, name: str, pop: str = None) -> Parameter:
 
         raise KeyError(f'Parameter "{name}" not found')
 
-    def sample(self, constant=True):
+    def sample(self, constant: bool = True, rng_sampler=None):
         """
         Return a sampled copy of the ParameterSet
 
         :param constant: If True, time series will be perturbed by a single constant offset. If False,
                          a different perturbation will be applied to each time specific value independently.
+        :param rng_sampler: Optional random number generator that may have been seeded to generate consistent results
         :return: A new :class:`ParameterSet` with perturbed values
 
         """
-
         new = sc.dcp(self)
-        for par in new.all_pars():
-            par.sample(constant)
+
+        for i, par in enumerate(new.all_pars()):
+            par.sample(constant=constant, rng_sampler=rng_sampler)
         return new
 
     def make_constant(self, year: float):
diff --git a/atomica/plotting.py b/atomica/plotting.py
index b0705cfc..db3c1375 100644
--- a/atomica/plotting.py
+++ b/atomica/plotting.py
@@ -516,6 +516,7 @@ def placeholder_pop():
                             numerator = sum(aggregated_outputs[x][output_name] * weights[x] for x in pop_labels)  # Add together all the outputs
                             denominator = sum([weights[x] for x in pop_labels])
                             vals = np.divide(numerator, denominator, out=np.full(numerator.shape, np.nan, dtype=float), where=numerator != 0)
+
                         else:
                             raise Exception(f'Unknown population aggregation method "{aggregation}"')
 
@@ -1132,7 +1133,9 @@ def set_colors(self, colors=None, results="all", pops="all", outputs="all", over
         elif isinstance(colors, list):
             assert len(colors) == len(targets), "Number of colors must either be a string, or a list with as many elements as colors to set"
             colors = colors
-        elif colors.startswith("#") or colors not in [m for m in plt.cm.datad if not m.endswith("_r")]:
+        elif isinstance(colors, tuple) and colors[0].startswith("#"): # tuple with (color, opacity)
+            colors = [colors for _ in range(len(targets))]  # Apply color to all requested outputs
+        elif isinstance(colors, str) and (colors.startswith("#") or colors not in [m for m in plt.cm.datad if not m.endswith("_r")]):
             colors = [colors for _ in range(len(targets))]  # Apply color to all requested outputs
         else:
             color_norm = matplotlib_colors.Normalize(vmin=-1, vmax=len(targets))
@@ -1658,7 +1661,8 @@ def process_input_stacks(input_stacks, available_items):
     return figs
 
 
-def plot_series(plotdata, plot_type="line", axis=None, data=None, legend_mode=None, lw=None, n_cols: int = None) -> list:
+def plot_series(plotdata, plot_type="line", axis=None, data=None, legend_mode=None, lw=None, n_cols: int = None,
+                colors=None) -> list:
     """
     Produce a time series plot
 
@@ -1671,6 +1675,7 @@ def plot_series(plotdata, plot_type="line", axis=None, data=None, legend_mode=No
     :param lw: override the default line width
     :param n_cols: If None (default), separate figures will be created for each axis. If provided, axes will be tiled as subplots in a single figure
                    window with the requested number of columns
+    :param colors: Colors to be passed to plotdata.set_colors
     :return: A list of newly created Figures
 
     """
@@ -1724,7 +1729,7 @@ def _prepare_figures(dim1, dim2, n_cols):
         logger.warning("At least one Series has only one timepoint. Series must have at least 2 time points to be rendered as a line - `plot_bars` may be more suitable for such data")
 
     if axis == "results":
-        plotdata.set_colors(results=plotdata.results.keys())
+        plotdata.set_colors(colors=colors, results=plotdata.results.keys())
 
         figs, axes = _prepare_figures(plotdata.pops, plotdata.outputs, n_cols)
 
@@ -1759,7 +1764,7 @@ def _prepare_figures(dim1, dim2, n_cols):
                 _render_legend(ax, plot_type)
 
     elif axis == "pops":
-        plotdata.set_colors(pops=plotdata.pops.keys())
+        plotdata.set_colors(colors=colors, pops=plotdata.pops.keys())
 
         figs, axes = _prepare_figures(plotdata.results, plotdata.outputs, n_cols)
 
@@ -1792,7 +1797,7 @@ def _prepare_figures(dim1, dim2, n_cols):
                 _render_legend(ax, plot_type)
 
     elif axis == "outputs":
-        plotdata.set_colors(outputs=plotdata.outputs.keys())
+        plotdata.set_colors(colors=colors, outputs=plotdata.outputs.keys())
 
         figs, axes = _prepare_figures(plotdata.results, plotdata.pops, n_cols)
 
diff --git a/atomica/programs.py b/atomica/programs.py
index 9b80b32c..2f4dd5ba 100644
--- a/atomica/programs.py
+++ b/atomica/programs.py
@@ -1141,7 +1141,7 @@ def get_outcomes(self, prop_coverage: dict) -> dict:
 
         return {(covout.par, covout.pop): covout.get_outcome(prop_coverage) for covout in self.covouts.values()}
 
-    def sample(self, constant: bool = True):
+    def sample(self, constant: bool = True, rng_sampler=None):
         """
         Perturb programs based on uncertainties
 
@@ -1158,9 +1158,9 @@ def sample(self, constant: bool = True):
 
         new = sc.dcp(self)
         for prog in new.programs.values():
-            prog.sample(constant)
+            prog.sample(constant, rng_sampler=rng_sampler)
         for covout in new.covouts.values():
-            covout.sample()
+            covout.sample(rng_sampler=rng_sampler)
         return new
 
 
@@ -1221,7 +1221,7 @@ def is_one_off(self) -> bool:
         """
         return "/year" not in self.unit_cost.units
 
-    def sample(self, constant: bool) -> None:
+    def sample(self, constant: bool, rng_sampler=None) -> None:
         """
         Perturb program values based on uncertainties
 
@@ -1231,11 +1231,12 @@ def sample(self, constant: bool) -> None:
 
         :param constant: If True, time series will be perturbed by a single constant offset. If False,
                          an different perturbation will be applied to each time specific value independently.
+        :param rng_sampler: Optional random number generator that may have been seeded to generate consistent results
         """
 
         self.spend_data = self.spend_data.sample(constant)
         self.unit_cost = self.unit_cost.sample(constant)
-        self.capacity_constraint = self.capacity_constraint.sample(constant)
+        self.capacity_constraint = self.capacity_constraint.sample(constant, rng_sampler)
         self.saturation = self.saturation.sample(constant)
         self.coverage = self.coverage.sample(constant)
 
@@ -1397,25 +1398,28 @@ def n_progs(self) -> int:
 
         return len(self.progs)
 
-    def sample(self) -> None:
+    def sample(self, rng_sampler=None) -> None:
         """
         Perturb the values entered in the databook
 
         The :class:`Covout` instance is modified in-place. Note that the program outcomes are scalars
         that do not vary over time - therefore, :meth:`Covout.sample()` does not have a ``constant``
         argument.
-
+        :param rng_sampler: Optional random number generator that may have been seeded to generate consistent results
         """
 
         if self.sigma is None:
             return
 
+        if rng_sampler is None:
+            rng_sampler = np.random
+
         for k, v in self.progs.items():
-            self.progs[k] = v + self.sigma * np.random.randn(1)[0]
+            self.progs[k] = v + self.sigma * rng_sampler.standard_normal(1)[0]
         # Perturb the interactions
         if self._interactions:
             for k, v in self.interactions.items():
-                self.interactions[k] = v + self.sigma * np.random.randn(1)[0]
+                self.interactions[k] = v + self.sigma * rng_sampler.standard_normal(1)[0]
             tokens = ["%s=%.4f" % ("+".join(k), v) for k, v in self.interactions.items()]
             self.imp_interaction = ",".join(tokens)
 
diff --git a/atomica/project.py b/atomica/project.py
index a7911dd4..3db746ac 100644
--- a/atomica/project.py
+++ b/atomica/project.py
@@ -44,10 +44,24 @@
 
 
 class ProjectSettings:
-    def __init__(self, sim_start=2000, sim_end=2035, sim_dt=0.25):
+    def __init__(self, sim_start:float=2000, sim_end:float=2035, sim_dt:float=0.25, multinomial:bool=False):
+        """
+        Project settings object
+
+        This contains various settings that are used by the Project when running simulations.
+        Most importantly, it contains the simulation start year, end year, and timestep, with
+        validation performed when setting these values.
+
+        :param sim_start: Simulation start year
+        :param sim_end: Simulation end year
+        :param sim_dt: Simulation time step
+        :param multinomial: Use multinomial mode with stochastic transitions
+        """
+
         self._sim_start = sim_start
         self._sim_dt = sim_dt
         self._sim_end = 0.0
+        self.multinomial = multinomial
         self.update_time_vector(end=sim_end)
 
     def __repr__(self):
@@ -93,9 +107,9 @@ def tvec(self) -> np.ndarray:
         :return: Array of simulation times
 
         """
-
         return np.arange(round((self.sim_end - self.sim_start) / self.sim_dt) + 1) * self.sim_dt + self.sim_start
 
+
     def update_time_vector(self, start: float = None, end: float = None, dt: float = None) -> None:
         """
         Update the project simulation times
@@ -474,7 +488,7 @@ def update_settings(self, sim_start=None, sim_end=None, sim_dt=None):
         """Modify the project settings, e.g. the simulation time vector."""
         self.settings.update_time_vector(start=sim_start, end=sim_end, dt=sim_dt)
 
-    def run_sim(self, parset=None, progset=None, progset_instructions=None, store_results=False, result_name: str = None):
+    def run_sim(self, parset=None, progset=None, progset_instructions=None, store_results=False, result_name: str = None, rng=None, acceptance_criteria=[]):
         """
         Run a single simulation
 
@@ -487,9 +501,14 @@ def run_sim(self, parset=None, progset=None, progset_instructions=None, store_re
         :param progset_instructions: A :class:`ProgramInstructions` instance. Programs will only be used if a instructions are provided
         :param store_results: If True, then the result will automatically be stored in ``self.results``
         :param result_name: Optionally assign a specific name to the result (otherwise, a unique default name will automatically be selected)
+        :param rng: Optionally a random number generator that may have been seeded to generate consistent results, or a random_seed used to generate a Generator
+        :param acceptance_criteria: Optional criteria to assert that outputs are within a given tolerance of a given value for given parameters over a given period of time
+        tolerance e.g.  = [{'parameter': 'incidence', 'population': 'adults', 'value': (0.02, 0.025), 't_range': (2018, 2018.99)}]
+
         :return: A :class:`Result` instance
 
         """
+        rng_sampler = rng if isinstance(rng, np.random._generator.Generator) else np.random.default_rng(rng)
 
         parset = self.parset(parset)
         if progset is not None:
@@ -508,14 +527,14 @@ def run_sim(self, parset=None, progset=None, progset_instructions=None, store_re
                 k += 1
 
         tm = sc.tic()
-        result = run_model(settings=self.settings, framework=self.framework, parset=parset, progset=progset, program_instructions=progset_instructions, name=result_name)
+        result = run_model(settings=self.settings, framework=self.framework, parset=parset, progset=progset, program_instructions=progset_instructions, name=result_name, rng_sampler=rng_sampler, acceptance_criteria=acceptance_criteria)
         logger.info('Elapsed time for running "%s": %ss', self.name, sc.sigfig(sc.toc(tm, output=True), 3))
         if store_results:
             self.results.append(result)
 
         return result
 
-    def run_sampled_sims(self, parset, progset=None, progset_instructions=None, result_names=None, n_samples: int = 1, parallel=False, max_attempts=None, num_workers=None) -> list:
+    def run_sampled_sims(self, parset, progset=None, progset_instructions=None, result_names=None, n_samples: int = 1, rand_seed: int = None, parallel=False, max_attempts=None, num_workers=None, acceptance_criteria=[]) -> list:
         """
         Run sampled simulations
 
@@ -539,6 +558,11 @@ def run_sampled_sims(self, parset, progset=None, progset_instructions=None, resu
         :param parallel: If True, run simulations in parallel (on Windows, must have ``if __name__ == '__main__'`` gating the calling code)
         :param max_attempts: Number of retry attempts for bad initializations
         :param num_workers: If ``parallel`` is True, this determines the number of parallel workers to use (default is usually number of CPUs)
+        :param rng_sampler: Optional random number generator that may have been seeded to generate consistent results
+        :param acceptance_criteria: Optional criteria to assert that outputs are within a given tolerance of a given value for given parameters in given years
+            tolerance e.g.  = [{'parameter': 'incidence', 'population': 'adults', 'value': (0.02, 0.025), 't': 2018}]
+
+
         :return: A list of Results that can be passed to `Ensemble.update()`. If multiple instructions are provided, the return value of this
                  function will be a list of lists, where the inner list iterates over different instructions for the same parset/progset samples.
                  It is expected in that case that the Ensemble's mapping function would take in a list of results
@@ -562,17 +586,27 @@ def run_sampled_sims(self, parset, progset=None, progset_instructions=None, resu
 
         show_progress = n_samples > 1 and logger.getEffectiveLevel() <= logging.INFO
 
+        if rand_seed is not None:
+            rng = np.random.default_rng(seed=rand_seed)
+            seed_samples = rng.integers(1e15, size=n_samples)
+        else:
+            seed_samples = [None] * n_samples
+
+        model_rngs = [np.random.default_rng(seed=seed) for seed in seed_samples]  # generate a RNG for each model
+
         if parallel:
-            fcn = functools.partial(_run_sampled_sim, proj=self, parset=parset, progset=progset, progset_instructions=progset_instructions, result_names=result_names, max_attempts=max_attempts)
-            results = parallel_progress(fcn, n_samples, show_progress=show_progress, num_workers=num_workers)
+            fcn = functools.partial(_run_sampled_sim, proj=self, parset=parset, progset=progset, progset_instructions=progset_instructions, result_names=result_names, max_attempts=max_attempts, acceptance_criteria=acceptance_criteria)
+            # as multiprocessing does not handle partial functions as compiled functions, need to send the rngs as kwargs in a dictionary, not as args to the partial function
+            model_rng_kwargs = [{"rng_sampler": rng} for rng in model_rngs]
+            results = parallel_progress(fcn, model_rng_kwargs, show_progress=show_progress, num_workers=num_workers)
         elif show_progress:
             # Print the progress bar if the logging level was INFO or lower
             # This means that the user can still set the logging level higher e.g. WARNING to suppress output from Atomica in general
             # (including any progress bars)
             with Quiet():
-                results = [_run_sampled_sim(self, parset, progset, progset_instructions, result_names, max_attempts=max_attempts) for _ in tqdm.trange(n_samples)]
+                results = [_run_sampled_sim(self, parset, progset, progset_instructions, result_names, max_attempts=max_attempts, rng_sampler=rng, acceptance_criteria=acceptance_criteria) for rng in tqdm.tqdm(model_rngs)]
         else:
-            results = [_run_sampled_sim(self, parset, progset, progset_instructions, result_names, max_attempts=max_attempts) for _ in range(n_samples)]
+            results = [_run_sampled_sim(self, parset, progset, progset_instructions, result_names, max_attempts=max_attempts, rng_sampler=rng, acceptance_criteria=acceptance_criteria) for rng in model_rngs]
 
         return results
 
@@ -745,7 +779,7 @@ def __setstate__(self, d):
         self.__dict__ = P.__dict__
 
 
-def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_names: list, max_attempts: int = None):
+def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_names: list, max_attempts: int = None, rng_sampler=None, acceptance_criteria=[]):
     """
     Internal function to run simulation with sampling
 
@@ -766,6 +800,9 @@ def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_n
     :param progset_instructions: A list of instructions to run against a single sample
     :param result_names: A list of result names (strings)
     :param max_attempts: Maximum number of sampling attempts before raising an error
+    :param rng_sampler: Optional random number generator that may have been seeded to generate consistent results
+    :param acceptance_criteria: Optional criteria to assert that outputs are within a given tolerance of a given value for given parameters over a given period of time
+        tolerance e.g.  = [{'parameter': 'incidence', 'population': 'adults', 'value': (0.02, 0.025), 't_range': (2018, 2018.99)}]
     :return: A list of results that either contains 1 result, or the same number of results as instructions
 
     """
@@ -775,16 +812,23 @@ def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_n
     if max_attempts is None:
         max_attempts = 50
 
+    if rng_sampler is None:
+        rng_sampler = np.random.default_rng()
+
+    # Set up separate seeds for each attempt and seed a rng for each one (in case scenarios diverge in later timepoints, it's important to START from the same seeds)
+    rand_seeds = rng_sampler.integers(1e15, size=max_attempts)
+
     attempts = 0
     while attempts < max_attempts:
         try:
+            rng_attempt = np.random.default_rng(seed=rand_seeds[attempts])
+            sampled_parset = parset.sample(rng_sampler=rng_attempt)
             if progset:
-                sampled_parset = parset.sample()
-                sampled_progset = progset.sample()
-                results = [proj.run_sim(parset=sampled_parset, progset=sampled_progset, progset_instructions=x, result_name=y) for x, y in zip(progset_instructions, result_names)]
+                sampled_progset = progset.sample(rng_sampler=rng_attempt)
+                results = [proj.run_sim(parset=sampled_parset, progset=sampled_progset, progset_instructions=x, result_name=y, rng=rng_attempt, acceptance_criteria=acceptance_criteria) for x, y in zip(progset_instructions, result_names)]
             else:
-                sampled_parset = parset.sample()
-                results = [proj.run_sim(parset=sampled_parset, result_name=y) for y in result_names]
+                results = [proj.run_sim(parset=sampled_parset, result_name=y, rng=rng_attempt, acceptance_criteria=acceptance_criteria) for y in result_names]
+
             return results
         except BadInitialization:
             attempts += 1
diff --git a/atomica/results.py b/atomica/results.py
index 9fe1ab3e..f3fa1df7 100644
--- a/atomica/results.py
+++ b/atomica/results.py
@@ -840,30 +840,21 @@ def _write_df(writer, formats, sheet_name, df, level_ordering):
     for level in level_ordering:
         order[level] = list(dict.fromkeys(df.index.get_level_values(level)))
 
-    required_width = [0] * (len(level_ordering) - 1)
-
-    row = 0
+    df = df.reorder_levels(level_ordering)
+    for i in range(len(level_ordering)):
+        df = df.reindex(order[level_ordering[i]], level=i)
 
     worksheet = writer.book.add_worksheet(sheet_name)
     writer.sheets[sheet_name] = worksheet  # Need to add it to the ExcelWriter for it to behave properly
 
-    for title, table in df.groupby(level=level_ordering[0], sort=False):
-        worksheet.write_string(row, 0, title, formats["center_bold"])
-        row += 1
-
-        table.reset_index(level=level_ordering[0], drop=True, inplace=True)  # Drop the title column
-        table = table.reorder_levels(level_ordering[1:])
-        for i in range(1, len(level_ordering)):
-            table = table.reindex(order[level_ordering[i]], level=i - 1)
-        table.index = table.index.set_names([level_substitutions[x] if x in level_substitutions else x.title() for x in table.index.names])
-        table.to_excel(writer, sheet_name=sheet_name, startcol=0, startrow=row, merge_cells=False)
-        row += table.shape[0] + 2
+    df.index = df.index.set_names([level_substitutions[x] if x in level_substitutions else x.title() for x in df.index.names])
+    df.to_excel(writer, sheet_name, startcol=0, startrow=0, merge_cells=False)
 
-        required_width[0] = max(required_width[0], len(title))
-        for i in range(0, len(required_width)):
-            required_width[i] = max(required_width[i], max(table.index.get_level_values(i).str.len()))
+    required_width = [len(name) for name in df.index.names]
+    for i in range(len(required_width)):
+        required_width[i] = max(len(val) for val in order[level_ordering[i]])
 
-    for i in range(0, len(required_width)):
+    for i in range(len(required_width)):
         if required_width[i] > 0:
             worksheet.set_column(i, i, required_width[i] * 1.1 + 1)
 
@@ -1231,7 +1222,7 @@ def plot_series(self, fig=None, style="quartile", results=None, outputs=None, po
 
         if fig is None:
             fig = plt.figure()
-        ax = plt.gca()
+        ax = fig.gca()
 
         series_lookup = self._get_series()
 
@@ -1266,7 +1257,7 @@ def plot_series(self, fig=None, style="quartile", results=None, outputs=None, po
 
             proposed_label = "%s (%s)" % (output, these_series[0].unit_string)
             if ax.yaxis.get_label().get_text():
-                assert proposed_label == ax.yaxis.get_label().get_text(), "The outputs being superimposed have different units"
+                assert proposed_label == ax.yaxis.get_label().get_text(), f"The outputs being superimposed have different units: {proposed_label} vs {ax.yaxis.get_label().get_text()}"
             else:
                 ax.set_ylabel(proposed_label)
 
diff --git a/atomica/system.py b/atomica/system.py
index 9e3c149c..ab11dd17 100644
--- a/atomica/system.py
+++ b/atomica/system.py
@@ -85,7 +85,7 @@ class FrameworkSettings:
     DEFAULT_SYMBOL_INAPPLICABLE = "N.A."
     DEFAULT_POP_TYPE = "default"
 
-    RESERVED_KEYWORDS = ["t", "flow", "all", "dt", "total"]  # A code_name in the framework cannot be equal to one of these values
+    RESERVED_KEYWORDS = ["t", "flow", "all", "dt", "total", "constant"]  # A code_name in the framework cannot be equal to one of these values
     RESERVED_KEYWORDS += supported_functions.keys()
 
     RESERVED_SYMBOLS = set(":,;/+-*'\" @")  # A code_name in the framework (for characs, comps, pars) cannot contain any of these characters
diff --git a/atomica/utils.py b/atomica/utils.py
index c20dfd7c..86bb47a6 100644
--- a/atomica/utils.py
+++ b/atomica/utils.py
@@ -595,7 +595,7 @@ def interpolate(self, t2: np.array, method="linear", **kwargs) -> np.array:
         interpolator = method(t1, v1, **kwargs)
         return interpolator(t2)
 
-    def sample(self, constant=True):
+    def sample(self, constant: bool = True, rng_sampler=None):
         """
         Return a sampled copy of the TimeSeries
 
@@ -604,6 +604,8 @@ def sample(self, constant=True):
 
         :param constant: If True, time series will be perturbed by a single constant offset. If False,
                          an different perturbation will be applied to each time specific value independently.
+        :param rng_sampler: Optional random number generator that may have been seeded to generate consistent results
+
         :return: A copied ``TimeSeries`` with perturbed values
 
         """
@@ -611,9 +613,12 @@ def sample(self, constant=True):
         if self._sampled:
             raise Exception("Sampling has already been performed - can only sample once")
 
+        if rng_sampler is None:
+            rng_sampler = np.random  #
+
         new = self.copy()
         if self.sigma is not None:
-            delta = self.sigma * np.random.randn(1)[0]
+            delta = self.sigma * rng_sampler.standard_normal(1)[0]
             if self.assumption is not None:
                 new.assumption += delta
 
@@ -622,7 +627,7 @@ def sample(self, constant=True):
                 new.vals = [v + delta for v in new.vals]
             else:
                 # Sample again for each data point
-                for i, (v, delta) in enumerate(zip(new.vals, self.sigma * np.random.randn(len(new.vals)))):
+                for i, (v, delta) in enumerate(zip(new.vals, self.sigma * rng_sampler.standard_normal(len(new.vals)))):
                     new.vals[i] = v + delta
 
         # Sampling flag only needs to be set if the TimeSeries had data to change
@@ -865,7 +870,10 @@ def callback(result, idx):
             jobs.append(pool.apply_async(fcn, callback=partial(callback, idx=i)))
     else:
         for i, x in enumerate(inputs):
-            jobs.append(pool.apply_async(fcn, args=(x,), callback=partial(callback, idx=i)))
+            if isinstance(x, dict):
+                pool.apply_async(fcn, args=(), kwds=x, callback=partial(callback, idx=i))  # if the list of inputs is given as a dictionary then pass as kwargs (kwds) instead
+            else:
+                pool.apply_async(fcn, args=(x,), callback=partial(callback, idx=i))
 
     pool.close()
     pool.join()
@@ -1008,6 +1016,20 @@ def stop_logging() -> None:
             # (not supposed to happen though) then we would want to close them all
 
 
+def stochastic_rounding(x, rng_sampler=None):
+    """
+    Stochastically round a float up or down to an integer-equivalent value (note: returns still as a float)
+    :param x: value to be rounded
+    """
+    sample = np.random.random() if rng_sampler is None else rng_sampler.random()
+
+    floor = np.floor(x)
+    remainder = x - floor
+    if remainder:
+        x = floor + int(sample < remainder)
+    return x
+
+
 def get_sigfigs_necessary(x, y, min_sigfigs: int = 2) -> int:
     """
     Get how many significant figures are necessary to tell the difference between two numbers
diff --git a/docs/examples/Uncertainty.ipynb b/docs/examples/Uncertainty.ipynb
index daef6d9d..fa652784 100644
--- a/docs/examples/Uncertainty.ipynb
+++ b/docs/examples/Uncertainty.ipynb
@@ -95,7 +95,7 @@
     "import matplotlib.pyplot as plt\n",
     "from scipy import stats\n",
     "import sciris as sc\n",
-    "np.random.seed(3)\n",
+    "np.random.seed(5)\n",
     "testdir = '../../tests/'\n",
     "\n",
     "P = at.Project(framework=testdir + 'test_uncertainty_framework.xlsx', databook=testdir + 'test_uncertainty_databook.xlsx')\n",
@@ -146,7 +146,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "If there is uncertainty, we can obtain a sampled epidemiological prediction by sampling from all of the quantities containing uncertinty. For example, replacing the value of `214` above with a sample from the distribution $\\ \\mathcal{N}(214,20)$ - along with all of the other quantities. The input to the model is effectively a parameter set derived from the original set of parameters, but with values replaced by sampled values where appropriate. You can obtain a sampled parameter set simply by calling the `ParameterSet.sample()` method:"
+    "If there is uncertainty, we can obtain a sampled epidemiological prediction by sampling from all of the quantities containing uncertainty. For example, replacing the value of `214` above with a sample from the distribution $\\ \\mathcal{N}(214,20)$ - along with all of the other quantities. The input to the model is effectively a parameter set derived from the original set of parameters, but with values replaced by sampled values where appropriate. You can obtain a sampled parameter set simply by calling the `ParameterSet.sample()` method:"
    ]
   },
   {
@@ -1035,9 +1035,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python [conda env:atomica37]",
+   "display_name": "Python [conda env:atomica312]",
    "language": "python",
-   "name": "conda-env-atomica37-py"
+   "name": "conda-env-atomica312-py"
   },
   "language_info": {
    "codemirror_mode": {
@@ -1049,7 +1049,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.3"
+   "version": "3.12.2"
   },
   "toc": {
    "base_numbering": 1,
diff --git a/tests/conftest.py b/tests/conftest.py
index 9934aed3..3439bfe3 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -6,3 +6,12 @@
 def close_figures():
     yield
     plt.close("all")
+
+
+# conftest.py
+def pytest_collection_modifyitems(items):
+    for i, item in enumerate(items):
+        if item.name == "test_model[tb]":
+            items.pop(i)
+            items.insert(0, item)
+            break
diff --git a/tests/test_tox_library.py b/tests/test_tox_library.py
index 110dd246..a5e59c66 100644
--- a/tests/test_tox_library.py
+++ b/tests/test_tox_library.py
@@ -173,7 +173,7 @@ def run_regenerated_framework(proj):
 # each function call as a separate test, which can be run in parallel
 # simply by adding `pytest -n 4` for instance. Or via tox
 # `tox -- -o -n 4` (because the default config is to use n=2 for TravisCI)
-@pytest.mark.parametrize("model", models)
+@pytest.mark.parametrize("model", at.demo_projects)
 def test_model(model):
 
     framework_file = at.LIBRARY_PATH / f"{model}_framework.xlsx"
diff --git a/tests/test_tox_markovchain.py b/tests/test_tox_markovchain.py
new file mode 100644
index 00000000..c526cea2
--- /dev/null
+++ b/tests/test_tox_markovchain.py
@@ -0,0 +1,121 @@
+import atomica as at
+import numpy as np
+import os
+import matplotlib.pyplot as plt
+import sciris as sc
+import pytest
+
+
+dirname = at.parent_dir()
+# np.seterr(all='raise')
+
+
+@pytest.mark.parametrize("demo", at.demo_projects)
+def test_functional_all_framework_databook(demo):
+    # Test that all demo frameworks can be run with stochastic mode turned on
+    P = at.demo(demo, do_run=False)
+    P.settings.multinomial = True
+    P.run_sim(result_name="DTMC")
+
+
+@pytest.mark.parametrize("demo", ["sir", "udt", "tb_simple"])  # This test is expensive so only run it for a few lightweight models
+def test_compare_deterministic_markovchain(demo):
+
+    seed = 0
+    N = 1000
+
+    y_cutoff = 10  # (deterministic) Compartment values below this will not be considered for the comparisons
+    index_cutoff = 1500
+    y_diff_cutoff = 0
+
+    ratio_cutoff = 0.05
+
+    def filter_points(t, det, mc, mc_mean):
+        # t_fil, det_fil, mc_fil = sc.dcp(t), sc.dcp(det), sc.dcp(mc)
+
+        indices = np.ones(t.shape)
+        indices = np.logical_and(indices, det >= y_cutoff)
+        indices = np.logical_and(indices, np.arange(len(indices)) <= index_cutoff)
+        indices = np.logical_and(indices, np.abs(det - mc_mean) > y_diff_cutoff)
+
+        return t[indices], det[indices], mc[:, indices], mc_mean[indices]
+
+    def get_array_outputs(results, output, pop):
+        extracted = [None] * len(results)
+        for i, res in enumerate(results):
+            extracted[i] = res.model.get_pop(pop).get_variable(output)[0].vals
+
+        array = np.stack(extracted, axis=0)
+        t = results[0].model.get_pop(pop).get_variable(output)[0].t
+        return t, array
+
+    def get_average_output(results, output, pop):
+        t, array = get_array_outputs(results, output, pop)
+        average = np.mean(array, axis=0)
+
+        return t, average, array
+
+    def do_hypothesis_test(det_array, mc_results):
+        import scipy
+
+        res = scipy.stats.ttest_1samp(mc_results, det_array, keepdims=True)
+        return res
+
+    def run_simple_percent_test(P, baseline_results, mc_results):
+        baseline_results = sc.promotetolist(baseline_results)
+        mc_results = sc.promotetolist(mc_results)
+
+        compartment_names = [name for name in baseline_results[0].model.framework.comps.index.values]
+
+        max_max_diff = 0
+
+        for compartment_name in compartment_names:
+            for pop in baseline_results[0].pop_names:
+                try:
+                    baseline_results[0].model.get_pop(pop).get_variable(compartment_name)[0].vals
+                except:
+                    continue  # This population doesn't have this compartment
+
+                t, mc_mean, mc_array = get_average_output(mc_results, compartment_name, pop)
+                _, det_mean, _ = get_average_output(baseline_results, compartment_name, pop)
+
+                t, det_mean, mc_array, mc_mean = filter_points(t, det_mean, mc_array, mc_mean)
+
+                if len(t) == 0:  # No filtered outputs so skip this output, pop combo
+                    continue
+
+                indices = det_mean != 0
+                ratio = np.ones(mc_mean.shape)
+                ratio[indices] = mc_mean[indices] / det_mean[indices]
+                if min(ratio) < 1 - ratio_cutoff or max(ratio) > 1 + ratio_cutoff:
+                    print("RATIO", compartment_name, pop, min(ratio), max(ratio), list(zip(ratio, det_mean, mc_mean)))
+                    indices = np.logical_or(ratio < 1 - ratio_cutoff, ration > 1 + ratio_cutoff)
+
+                    raise Exception(f"The Markov Chain runs (N={len(mc_results)} differed too much from the baseline runs (N={len(baseline_results)} at starting at timestep={indices[0]}, t={t[indices[0]]}, with compartment={compartment_name}, population={pop}. Baseline={det_mean[indices[0]]}, MC={mc_mean[indices[0]]}, difference={100*(mc_mean[indices[0]]/det_mean[indices[0]] - 1):.2f}")
+
+                max_diff = np.max(np.abs(ratio - 1))
+                max_max_diff = max(max_diff, max_max_diff)
+
+        print(f"Success: max difference {max_max_diff*100:.2f}%")
+
+    seed = np.random.default_rng(seed).integers(2**32 - 1)
+
+    P = at.demo(demo, do_run=False)
+    res = P.run_sim(result_name="Original")
+
+    baseline_results = [res]
+
+    P.settings.multinomial = True
+    mc_results = [None] * N
+    with at.Quiet():  # Using at.Quiet should automatically reset the logging level if an exception occurs
+        for i in range(N):
+            this_seed = (i + seed + hash(demo)) % (2**32)
+            mc_results[i] = P.run_sim(result_name=f"DTMC {i}", rng=this_seed)
+
+    run_simple_percent_test(P, baseline_results=baseline_results, mc_results=mc_results)
+
+
+if __name__ == "__main__":
+    for demo in at.demo_projects:
+        test_functional_all_framework_databook()
+    test_compare_deterministic_markovchain("sir")
diff --git a/tests/test_uncertainty_framework.xlsx b/tests/test_uncertainty_framework.xlsx
index 120b1735..37073ea5 100644
Binary files a/tests/test_uncertainty_framework.xlsx and b/tests/test_uncertainty_framework.xlsx differ