diff --git a/examples/highlevel_models/02-highlevel-models.py b/examples/highlevel_models/02-highlevel-models.py
new file mode 100644
index 000000000..13cdf5f9d
--- /dev/null
+++ b/examples/highlevel_models/02-highlevel-models.py
@@ -0,0 +1,52 @@
+import stormpy
+
+import stormpy.examples
+import stormpy.examples.files
+
+
+def example_highlevel_models():
+    """
+    This method considers a grid with a constant V.
+    We create an all-in-one MDP, i.e., we create an MDP that contains copies for different values of V.
+    By model checking, we obtain values for different values of V.
+    """
+    path = stormpy.examples.files.prism_mdp_slipgrid_sketch
+    prism_program = stormpy.parse_prism_program(path)
+    prop_string = "Rmin=? [F \"target\"]"
+    properties = stormpy.parse_properties(prop_string, prism_program)
+    # Note: The following is necessary for the code to work correctly, as otherwise there will be a default initial value injected.
+    prism_program = prism_program.replace_variable_initialization_by_init_expression()
+
+    emgr = prism_program.expression_manager
+    Vvar = prism_program.get_constant("V")
+    Vvalues = {2,4,6,8,10,12}
+    all_in_one_program = prism_program.replace_constant_by_variable(Vvar, emgr.create_integer(min(Vvalues)), emgr.create_integer(max(Vvalues)))
+    options = [stormpy.Expression.Eq(Vvar.expression_variable.get_expression(), emgr.create_integer(val)) for val in Vvalues]
+    all_in_one_program.update_initial_states_expression(stormpy.Expression.And(all_in_one_program.initial_states_expression, stormpy.Expression.Disjunction(options)))
+    print(all_in_one_program)
+
+    for v in all_in_one_program.constants:
+        print(f"Variable {v.name} is {v.definition}")
+    build_options = stormpy.BuilderOptions([p.raw_formula for p in properties])
+    build_options.set_build_state_valuations()
+    sparse_model = stormpy.build_sparse_model_with_options(all_in_one_program, build_options)
+    print(sparse_model)
+    result = stormpy.check_model_sparse(sparse_model, properties[0], extract_scheduler=True, only_initial_states=True)
+    filter_sparse = stormpy.create_filter_initial_states_sparse(sparse_model)
+    result.filter(filter_sparse)
+    result_for_vvar = {int(sparse_model.state_valuations.get_integer_value(x, Vvar.expression_variable)) : result.at(x) for x in sparse_model.initial_states}
+    print(result_for_vvar)
+
+    dd_model = stormpy.build_symbolic_model(all_in_one_program)
+    print(dd_model)
+    filter_dd = stormpy.create_filter_initial_states_symbolic(dd_model)
+    result = stormpy.check_model_dd(dd_model, properties[0], only_initial_states=True)
+    result.filter(filter_dd)
+
+    resultvalues = result.get_values()
+    # Notice that the manager seems to be different than the manager for the prism program. So we solve it like this
+    result_for_vvar = {int(valuation.get_integer_value(valuation.expression_manager.get_variable("V"))) : mc_result for valuation, mc_result in resultvalues }
+    print(result_for_vvar)
+
+if __name__ == '__main__':
+    example_highlevel_models()
\ No newline at end of file
diff --git a/lib/stormpy/examples/files.py b/lib/stormpy/examples/files.py
index d7e8167fb..8b75d0247 100644
--- a/lib/stormpy/examples/files.py
+++ b/lib/stormpy/examples/files.py
@@ -38,6 +38,8 @@ def _path(folder, file):
 prism_mdp_maze = _path("mdp", "maze_2.nm")
 """Prism example for the maze MDP"""
 prism_mdp_slipgrid = _path("mdp", "slipgrid.nm")
+"""Prism example for the maze MDP towards sketching"""
+prism_mdp_slipgrid_sketch = _path("mdp", "slipgrid_sketch.nm")
 """Prism example for a slippery grid with fixed dimensions"""
 drn_pomdp_maze = _path("pomdp", "maze.drn")
 """DRN example for the maze POMDP"""
diff --git a/lib/stormpy/examples/files/mdp/slipgrid_sketch.nm b/lib/stormpy/examples/files/mdp/slipgrid_sketch.nm
new file mode 100644
index 000000000..1b1470fff
--- /dev/null
+++ b/lib/stormpy/examples/files/mdp/slipgrid_sketch.nm
@@ -0,0 +1,31 @@
+mdp
+
+const double p = 0.4;
+const int N = 10;
+const int M = 10;
+const int V;
+
+module grid
+    x : [1..N] init 1;
+    y : [1..M] init 1;
+    [north] x > 1 -> 1-p: (x'=x-1) + p: (x'=x);
+    [south] x < V -> 1-p: (x'=x+1) + p: (x'=x);
+    [east] y > 1 -> 1-p: (y'=y-1) + p: (x'=x);
+    [west] y < M -> 1-p: (y'=y+1) + p: (x'=x);
+endmodule
+
+// reward structure (number of steps to reach the target)
+rewards
+
+	[east] true : 1;
+	[west] true : 1;
+	[north] true : 1;
+	[south] true : V;
+
+endrewards
+
+// target observation
+label "target" = x=3 & y=3;
+label "goal" = x=N & y=M;
+
+
diff --git a/src/storage/prism.cpp b/src/storage/prism.cpp
index 076fedaf6..fe411f4fe 100644
--- a/src/storage/prism.cpp
+++ b/src/storage/prism.cpp
@@ -3,6 +3,7 @@
 #include <boost/variant.hpp>
 #include <boost/algorithm/string/join.hpp>
 #include <random>
+#include <optional>
 #include "src/helpers.h"
 #include <storm/storage/expressions/ExpressionManager.h>
 #include <storm/storage/expressions/ExpressionEvaluator.h>
@@ -48,6 +49,11 @@ void define_prism(py::module& m) {
             .def("substitute_nonstandard_predicates", &Program::substituteNonStandardPredicates, "Remove nonstandard predicates from the prism program")
             .def("used_constants",&Program::usedConstants, "Compute Used Constants")
             .def("label_unlabelled_commands", &Program::labelUnlabelledCommands, "Label unlabelled commands", py::arg("name_suggestions"))
+            .def("replace_variable_initialization_by_init_expression", &Program::replaceVariableInitializationByInitExpression, "Replaces initializations from the individual variables to a init expression.")
+            .def("replace_constant_by_variable", &Program::replaceConstantByVariable, py::arg("constant"), py::arg("lower_bound"), py::arg("upper_bound"), py::arg("observable") = true, "Operation to take a constant and make it into a global variable (ranging from lower to upper bound).")
+            .def_property_readonly("has_initial_states_expression", &Program::hasInitialConstruct, "Is an initial states expression given.")
+            .def_property_readonly("initial_states_expression", [](Program const& p) {return p.hasInitialConstruct() ? std::make_optional(p.getInitialStatesExpression()) : std::nullopt;}, "Get the initial states expression, or none, if none exists")
+            .def("update_initial_states_expression", &Program::updateInitialStatesExpression, py::arg("new_expression"), "Replace initial expression. Can only be called if initial expression exists.")
             .def("has_constant", &Program::hasConstant, py::arg("name"))
             .def("get_constant", &Program::getConstant, py::arg("name"), "Requires that the program has a constant with this name")
             .def("has_reward_model", [](Program const& p, std::string const& name) {return p.hasRewardModel(name);},  py::arg("name"), "Is a reward model with the specified name defined?")