diff --git a/python/sdist/amici/debugging/__init__.py b/python/sdist/amici/debugging/__init__.py
new file mode 100644
index 0000000000..81663d17b9
--- /dev/null
+++ b/python/sdist/amici/debugging/__init__.py
@@ -0,0 +1,45 @@
+"""Functions for debugging AMICI simulation failures."""
+import amici
+import numpy as np
+
+
+def get_model_for_preeq(model: amici.Model, edata: amici.ExpData):
+    """Get a model set-up to simulate the preequilibration condition as
+    specified in `edata`.
+
+    Useful for analyzing simulation errors during preequilibration.
+    Simulating the returned model will reproduce the behavior of
+    simulation-based preequilibration.
+
+    Note that for models with events, the simulation results may differ.
+    During preequilibration, event-handling is disabled. However, when
+    simulating the returned model, event-handling will be enabled.
+    For events triggered at fixed timepoints, this can be avoided by setting
+    :meth:`t0 <amici.Model.setT0>` to a timepoints after the last trigger
+    timepoint.
+
+    :param model:
+        The model for which *edata* was generated.
+    :param edata:
+        The experimental data object with a preequilibration condition
+        specified.
+    :return:
+        A copy of *model* with the same parameters, initial states, ...
+        as *amici_model* for the preequilibration condition.
+        Output timepoints are set to ``[inf]`` and will have to be adjusted.
+    """
+    model = model.clone()
+    model.setTimepoints([np.inf])
+    model.setFixedParameters(edata.fixedParametersPreequilibration)
+    if edata.pscale:
+        model.setParameterScale(edata.pscale)
+    if edata.parameters:
+        model.setParameters(edata.parameters)
+    if edata.plist:
+        model.setParameterList(edata.plist)
+    model.setInitialStates(edata.x0)
+    # has to be set *after* parameter list/scale!
+    model.setInitialStateSensitivities(edata.sx0)
+    model.setT0(edata.tstart_)
+
+    return model
diff --git a/python/tests/test_preequilibration.py b/python/tests/test_preequilibration.py
index a42bc6354d..be447b0c54 100644
--- a/python/tests/test_preequilibration.py
+++ b/python/tests/test_preequilibration.py
@@ -5,7 +5,8 @@
 import amici
 import numpy as np
 import pytest
-from numpy.testing import assert_allclose
+from amici.debugging import get_model_for_preeq
+from numpy.testing import assert_allclose, assert_equal
 from test_pysb import get_data
 
 
@@ -633,3 +634,31 @@ def test_simulation_errors(preeq_fixture):
         assert rdata._swigptr.messages[2].identifier == "OTHER"
         assert rdata._swigptr.messages[3].severity == amici.LogSeverity_debug
         assert rdata._swigptr.messages[3].identifier == "BACKTRACE"
+
+
+def test_get_model_for_preeq(preeq_fixture):
+    (
+        model,
+        solver,
+        edata,
+        edata_preeq,
+        edata_presim,
+        edata_sim,
+        pscales,
+        plists,
+    ) = preeq_fixture
+    model.setSteadyStateSensitivityMode(
+        amici.SteadyStateSensitivityMode.integrationOnly
+    )
+    model_preeq = get_model_for_preeq(model, edata)
+    # the exactly same settings are used, so results should match exactly
+    rdata1 = amici.runAmiciSimulation(model_preeq, solver)
+    rdata2 = amici.runAmiciSimulation(model, solver, edata_preeq)
+    assert_equal(
+        rdata1.x,
+        rdata2.x,
+    )
+    assert_equal(
+        rdata1.sx,
+        rdata2.sx,
+    )