diff --git a/src/core/pla.cpp b/src/core/pla.cpp
index f554cd8..641aa9c 100644
--- a/src/core/pla.cpp
+++ b/src/core/pla.cpp
@@ -12,6 +12,17 @@ storm::modelchecker::parametric::RegionCheckResult checkRegion(std::shared_ptr<P
     return checker->analyzeRegion(region, initialResult, sampleVertices);
 }
 
+std::set<storm::Polynomial> gatherDerivatives(storm::models::sparse::Dtmc<storm::RationalFunction> const& model, carl::Variable const& var) {
+    std::set<storm::Polynomial> derivatives;
+    for (auto it : model.getTransitionMatrix()) {
+        storm::Polynomial pol = it.getValue().derivative(var, 1).nominator();
+        if (!pol.isConstant()) {
+            derivatives.insert(pol);
+        }
+    }
+    return derivatives;
+}
+
 // Define python bindings
 void define_pla(py::module& m) {
     
@@ -44,4 +55,6 @@ void define_pla(py::module& m) {
         .def("specify_formula", &specifyFormula, "Specify formula", py::arg("formula"))
         .def("check_region", &checkRegion, "Check region", py::arg("region"), py::arg("initialResult") = storm::modelchecker::parametric::RegionCheckResult::Unknown, py::arg("sampleVertices") = false)
     ;
+
+    m.def("gather_derivatives", &gatherDerivatives, "Gather all derivatives of transition probabilities", py::arg("model"), py::arg("var"));
 }
diff --git a/tests/core/test_pla.py b/tests/core/test_pla.py
index be48f7a..4537d40 100644
--- a/tests/core/test_pla.py
+++ b/tests/core/test_pla.py
@@ -43,3 +43,18 @@ class TestModelChecking:
         region = stormpy.ParameterRegion("0.1<=pL<=0.73,0.2<=pK<=0.715", parameters)
         result = checker.check_region(region)
         assert result == stormpy.RegionCheckResult.ALLVIOLATED
+    
+    def test_derivatives(self):
+        import pycarl
+        program = stormpy.parse_prism_program(get_example_path("pdtmc", "brp16_2.pm"))
+        prop = "P<=0.84 [F s=5 ]"
+        formulas = stormpy.parse_properties_for_prism_program(prop, program)
+        model = stormpy.build_parametric_model(program, formulas)
+        assert model.nr_states == 613
+        assert model.nr_transitions == 803
+        assert model.model_type == stormpy.ModelType.DTMC
+        assert model.has_parameters
+        parameters = model.collect_probability_parameters()
+        assert len(parameters) == 2
+        derivatives = stormpy.gather_derivatives(model, list(parameters)[0])
+        assert len(derivatives) == 0