diff --git a/examples/tutorials/config_sample_2pop_vaccine_scenarios.yml b/examples/tutorials/config_sample_2pop_vaccine_scenarios.yml
new file mode 100644
index 000000000..28e7e9eae
--- /dev/null
+++ b/examples/tutorials/config_sample_2pop_vaccine_scenarios.yml
@@ -0,0 +1,168 @@
+name: sample_2pop
+setup_name: minimal
+start_date: 2020-02-01
+end_date: 2020-08-31
+nslots: 10
+
+subpop_setup:
+  geodata: model_input/geodata_sample_2pop.csv
+  mobility: model_input/mobility_sample_2pop.csv
+
+initial_conditions:
+  method: SetInitialConditions
+  initial_conditions_file: model_input/ic_2pop.csv
+  allow_missing_subpops: TRUE
+  allow_missing_compartments: TRUE
+
+compartments:
+  infection_stage: ["S", "E", "I", "R", "V"]
+
+seir:
+  integration:
+    method: rk4
+    dt: 1
+  parameters:
+    sigma:
+      value: 1 / 4
+    gamma:
+      value: 1 / 5
+    Ro:
+      value: 
+        distribution: truncnorm
+        mean: 2.5
+        sd: 0.1
+        a: 1.1
+        b: 3
+    omega_pess:
+        value: 0.02
+    omega_opt:
+        value: 0.01
+    nu_pess: 
+        value: 0.01
+    nu_opt: 
+        value: 0.03
+  transitions:
+    #infections
+    - source: ["S"]
+      destination: ["E"]
+      rate: ["Ro * gamma"]
+      proportional_to: [["S"],["I"]]
+      proportion_exponent: ["1","1"]
+    # progression to infectiousness
+    - source: ["E"]
+      destination: ["I"]
+      rate: ["sigma"]
+      proportional_to: ["E"]
+      proportion_exponent: ["1"]
+    # recovery
+    - source: ["I"]
+      destination: ["R"]
+      rate: ["gamma"]
+      proportional_to: ["I"]
+      proportion_exponent: ["1"]
+    #vaccination (offers complete protection)
+    - source: ["S"]
+      destination: ["V"]
+      rate: ["nu_pess + nu_opt"]
+      proportional_to: ["S"]
+      proportion_exponent: ["1"]
+    # waning of vaccine-induced protection
+    - source: ["V"]
+      destination: ["S"]
+      rate: ["omega_pess + omega_opt"]
+      proportional_to: ["V"]
+      proportion_exponent: ["1"]
+
+seir_modifiers:
+  scenarios:
+    - no_vax
+    - pess_vax
+    - opt_vax
+  modifiers:
+    pess_vax_nu: # turn off nu_opt, only nu_pess left
+      method: SinglePeriodModifier
+      parameter: nu_opt
+      period_start_date: 2020-02-01
+      period_end_date: 2020-08-31
+      subpop: "all"
+      value: 0
+    pess_vax_wane: # turn off omega_opt, only omega_pess left
+      method: SinglePeriodModifier
+      parameter: omega_opt
+      period_start_date: 2020-02-01
+      period_end_date: 2020-08-31
+      subpop: "all"
+      value: 0
+    pess_vax: # turn off all vaccination
+      method: StackedModifier
+      modifiers: ["pess_vax_nu","pess_vax_wane"]
+    opt_vax_nu: # turn off nu_pess, only nu_opt left
+      method: SinglePeriodModifier
+      parameter: nu_pess
+      period_start_date: 2020-02-01
+      period_end_date: 2020-08-31
+      subpop: "all"
+      value: 0
+    opt_vax_wane: # turn off omega_pess, only omega_opt left
+      method: SinglePeriodModifier
+      parameter: omega_pess
+      period_start_date: 2020-02-01
+      period_end_date: 2020-08-31
+      subpop: "all"
+      value: 0
+    opt_vax: # turn off all vaccination
+      method: StackedModifier
+      modifiers: ["opt_vax_nu","opt_vax_wane"]
+    no_vax: # turn off all vaccination
+      method: StackedModifier
+      modifiers: ["pess_vax","opt_vax"]
+
+
+
+outcomes:
+  method: delayframe
+  outcomes:
+    incidCase: #incidence of detected cases
+      source:
+        incidence:
+          infection_stage: "I"
+      probability:
+        value: 
+          distribution: truncnorm
+          mean: 0.5
+          sd: 0.05
+          a: 0.3
+          b: 0.7
+      delay:
+        value: 5
+    incidHosp: #incidence of hospitalizations
+      source:
+        incidence:
+          infection_stage: "I"
+      probability:
+        value: 0.05
+      delay:
+        value: 7
+      duration:
+        value: 10
+        name: currHosp # will track number of current hospitalizations (ie prevalence)
+    incidDeath: #incidence of deaths
+      source: incidHosp
+      probability:
+        value: 0.2
+      delay:
+        value: 14
+
+# outcome_modifiers:
+#   scenarios:
+#     - test_limits
+#   modifiers:
+#     # assume that due to limitations in testing, initially the case detection probability was lower
+#     test_limits:
+#       method: SinglePeriodModifier
+#       parameter: incidCase::probability
+#       subpop: "all"
+#       period_start_date: 2020-02-01
+#       period_end_date: 2020-06-01
+#       value: 0.5
+
diff --git a/flepimop/gempyor_pkg/tests/simulate/test_random_seir_parameter_draw_per_slot.py b/flepimop/gempyor_pkg/tests/simulate/test_random_seir_parameter_draw_per_slot.py
new file mode 100644
index 000000000..c12235195
--- /dev/null
+++ b/flepimop/gempyor_pkg/tests/simulate/test_random_seir_parameter_draw_per_slot.py
@@ -0,0 +1,78 @@
+import os
+from pathlib import Path
+import shutil
+
+from click.testing import CliRunner
+import pandas as pd
+import pytest
+
+from gempyor.simulate import _click_simulate
+
+
+@pytest.fixture
+def setup_sample_2pop_vaccine_scenarios(tmp_path: Path) -> Path:
+    tutorials_path = Path(os.path.dirname(__file__) + "/../../../../examples/tutorials")
+    for file in (
+        "config_sample_2pop_vaccine_scenarios.yml",
+        "model_input/geodata_sample_2pop.csv",
+        "model_input/mobility_sample_2pop.csv",
+        "model_input/ic_2pop.csv",
+    ):
+        source = tutorials_path / file
+        destination = tmp_path / file
+        destination.parent.mkdir(parents=True, exist_ok=True)
+        shutil.copy(source, destination)
+    return tmp_path
+
+
+def test_random_seir_parameter_draw_per_slot(
+    monkeypatch: pytest.MonkeyPatch, setup_sample_2pop_vaccine_scenarios: Path
+) -> None:
+    # Test setup
+    monkeypatch.chdir(setup_sample_2pop_vaccine_scenarios)
+
+    # Test execution of `gempyor-simulate`
+    runner = CliRunner()
+    result = runner.invoke(_click_simulate, ["config_sample_2pop_vaccine_scenarios.yml"])
+    assert result.exit_code == 0
+
+    # Get the contents of 'spar' and 'hpar' directories as DataFrames
+    spar_directory: Path | None = None
+    hpar_directory: Path | None = None
+    for p in setup_sample_2pop_vaccine_scenarios.rglob("*"):
+        if p.is_dir() and p.name == "spar":
+            spar_directory = p
+        elif p.is_dir() and p.name == "hpar":
+            hpar_directory = p
+        if spar_directory is not None and hpar_directory is not None:
+            break
+
+    def read_directory(directory: Path) -> list[pd.DataFrame]:
+        dfs: list[pd.DataFrame] | pd.DataFrame = []
+        for i, f in enumerate(sorted(list(directory.glob("*.parquet")))):
+            dfs.append(pd.read_parquet(f))
+            dfs[-1]["slot"] = i
+        dfs = pd.concat(dfs)
+        return dfs
+
+    hpar = read_directory(hpar_directory)
+    spar = read_directory(spar_directory)
+
+    # Test contents of 'spar'/'hpar' DataFrames
+    assert (
+        hpar[
+            (hpar["subpop"] == "large_province")
+            & (hpar["quantity"] == "probability")
+            & (hpar["outcome"] == "incidCase")
+        ]["value"].nunique()
+        == 10
+    )
+    assert (
+        hpar[
+            (hpar["subpop"] == "small_province")
+            & (hpar["quantity"] == "probability")
+            & (hpar["outcome"] == "incidCase")
+        ]["value"].nunique()
+        == 10
+    )
+    assert spar[spar["parameter"] == "Ro"]["value"].nunique() == 10