diff --git a/lib/stormpy/__init__.py b/lib/stormpy/__init__.py
index f70697c..67e998f 100644
--- a/lib/stormpy/__init__.py
+++ b/lib/stormpy/__init__.py
@@ -108,9 +108,11 @@ def perform_bisimulation(model, properties, bisimulation_type):
 
 def model_checking(model, property):
     if model.supports_parameters:
-        return core._parametric_model_checking(model, property.raw_formula)
+        task = core.ParametricCheckTask(property.raw_formula, False)
+        return core._parametric_model_checking_sparse_engine(model, task)
     else:
-        return core._model_checking(model, property.raw_formula)
+        task = core.CheckTask(property.raw_formula, False)
+        return core._model_checking_sparse_engine(model, task)
 
 
 def compute_prob01_states(model, phi_states, psi_states):
diff --git a/src/core/analysis.cpp b/src/core/analysis.cpp
index 2dda203..0684969 100644
--- a/src/core/analysis.cpp
+++ b/src/core/analysis.cpp
@@ -10,5 +10,4 @@ void define_graph_constraints(py::module& m) {
             .def_property_readonly("graph_preserving_constraints", &storm::analysis::ConstraintCollector<storm::RationalFunction>::getGraphPreservingConstraints)
     ;
 
-
 }
diff --git a/src/core/analysis.h b/src/core/analysis.h
index 379acd9..289a7ef 100644
--- a/src/core/analysis.h
+++ b/src/core/analysis.h
@@ -1,5 +1,4 @@
 #pragma once
 #include "common.h"
 
-
-void define_graph_constraints(py::module& m);
\ No newline at end of file
+void define_graph_constraints(py::module& m);
diff --git a/src/core/bisimulation.cpp b/src/core/bisimulation.cpp
index cdd3bf0..7feddd0 100644
--- a/src/core/bisimulation.cpp
+++ b/src/core/bisimulation.cpp
@@ -1,19 +1,13 @@
 #include "bisimulation.h"
 
-// Thin wrapper for bisimulation
-template<typename ValueType>
-std::shared_ptr<storm::models::sparse::Model<ValueType>> performBisimulation(std::shared_ptr<storm::models::sparse::Model<ValueType>> const& model, std::vector<std::shared_ptr<storm::logic::Formula const>> const& formulas, storm::storage::BisimulationType bisimulationType) {
-    return storm::performBisimulationMinimization<storm::models::sparse::Model<ValueType>>(model, formulas, bisimulationType);
-}
-
 // Define python bindings
 void define_bisimulation(py::module& m) {
 
     // Bisimulation
-    m.def("_perform_bisimulation", &performBisimulation<double>, "Perform bisimulation", py::arg("model"), py::arg("formulas"), py::arg("bisimulation_type"));
-    m.def("_perform_parametric_bisimulation", &performBisimulation<storm::RationalFunction>, "Perform bisimulation on parametric model", py::arg("model"), py::arg("formulas"), py::arg("bisimulation_type"));
+    m.def("_perform_bisimulation", &storm::api::performBisimulationMinimization<double>, "Perform bisimulation", py::arg("model"), py::arg("formulas"), py::arg("bisimulation_type"));
+    m.def("_perform_parametric_bisimulation", &storm::api::performBisimulationMinimization<storm::RationalFunction>, "Perform bisimulation on parametric model", py::arg("model"), py::arg("formulas"), py::arg("bisimulation_type"));
 
-    // BisimulationType 
+    // BisimulationType
     py::enum_<storm::storage::BisimulationType>(m, "BisimulationType", "Types of bisimulation")
         .value("STRONG", storm::storage::BisimulationType::Strong)
         .value("WEAK", storm::storage::BisimulationType::Weak)
diff --git a/src/core/common.h b/src/core/common.h
index e8466b2..e961661 100644
--- a/src/core/common.h
+++ b/src/core/common.h
@@ -1,3 +1,3 @@
 #include "src/common.h"
-#include "storm/utility/storm.h"
+#include "storm/api/storm.h"
 #include <pybind11/stl.h>
diff --git a/src/core/core.cpp b/src/core/core.cpp
index 4a7dcf1..fa6fedc 100644
--- a/src/core/core.cpp
+++ b/src/core/core.cpp
@@ -1,5 +1,7 @@
 #include "core.h"
+#include "storm/utility/initialize.h"
 #include "storm/utility/DirectEncodingExporter.h"
+#include "storm/storage/ModelFormulasPair.h"
 
 void define_core(py::module& m) {
     // Init
@@ -12,7 +14,9 @@ void define_core(py::module& m) {
 
 void define_parse(py::module& m) {
     // Parse formulas
-    m.def("parse_properties", &storm::parsePropertiesForExplicit, R"dox(
+    m.def("parse_properties", [](std::string const& inputString, boost::optional<std::set<std::string>> const& propertyFilter = boost::none) {
+                return storm::api::parseProperties(inputString, propertyFilter);
+            }, R"dox(
           
           Parse properties given in the prism format.
           
@@ -20,7 +24,7 @@ void define_parse(py::module& m) {
           :param str property_filter: A filter
           :return: A list of properties
           )dox", py::arg("formula_string"), py::arg("property_filter") = boost::none);
-    m.def("parse_properties_for_prism_program", &storm::parsePropertiesForPrismProgram,  R"dox(
+    m.def("parse_properties_for_prism_program", &storm::api::parsePropertiesForPrismProgram,  R"dox(
           
           Parses properties given in the prism format, allows references to variables in the prism program.
           
@@ -39,28 +43,21 @@ void define_parse(py::module& m) {
                 return pair.formulas;
             }, "The formulas")
     ;
-
-    // Parse explicit models
-    m.def("parse_explicit_model", &storm::parser::AutoParser<>::parseModel, "Parse explicit model", py::arg("transition_file"), py::arg("labeling_file"), py::arg("state_reward_file") = "", py::arg("transition_reward_file") = "", py::arg("choice_labeling_file") = "");
 }
 
 // Thin wrapper for model building using one formula as argument
 template<typename ValueType>
-std::shared_ptr<storm::models::ModelBase> buildSymbolicModel(storm::prism::Program const& program, std::shared_ptr<storm::logic::Formula const> const& formula) {
-    return storm::buildSymbolicModel<ValueType>(program, std::vector<std::shared_ptr<storm::logic::Formula const>>(1,formula));
-}
-
-template<typename ValueType>
-std::shared_ptr<storm::models::ModelBase> buildSparseModel(storm::prism::Program const& program, std::shared_ptr<storm::logic::Formula const> const& formula) {
-    return storm::buildSparseModel<ValueType>(program, std::vector<std::shared_ptr<storm::logic::Formula const>>(1,formula));
+std::shared_ptr<storm::models::ModelBase> buildSparseModel(storm::storage::SymbolicModelDescription const& modelDescription, std::vector<std::shared_ptr<storm::logic::Formula const>> const& formulas, bool buildChoiceLabels = false, bool buildChoiceOrigins = false) {
+    return storm::api::buildSparseModel<ValueType>(modelDescription, formulas, buildChoiceLabels, buildChoiceOrigins);
 }
 
 void define_build(py::module& m) {
     // Build model
-    m.def("_build_sparse_model_from_prism_program", &storm::buildSparseModel<double>, "Build the model", py::arg("program"), py::arg("formulas") = std::vector<std::shared_ptr<storm::logic::Formula const>>());
-    m.def("_build_sparse_parametric_model_from_prism_program", &storm::buildSparseModel<storm::RationalFunction>, "Build the parametric model", py::arg("program"), py::arg("formulas"));
-    m.def("_build_sparse_model_from_drn", &storm::buildExplicitDRNModel<double>, "Build the model from DRN", py::arg("file"));
-    m.def("_build_sparse_parametric_model_from_drn", &storm::buildExplicitDRNModel<storm::RationalFunction>, "Build the parametric model from DRN", py::arg("file"));
+    m.def("_build_sparse_model_from_prism_program", &buildSparseModel<double>, "Build the model", py::arg("model_description"), py::arg("formulas") = std::vector<std::shared_ptr<storm::logic::Formula const>>(), py::arg("build_choice_labels") = false, py::arg("build_choice_origins") = false);
+    m.def("_build_sparse_parametric_model_from_prism_program", &buildSparseModel<storm::RationalFunction>, "Build the parametric model", py::arg("model_description"), py::arg("formulas"), py::arg("build_choice_labels") = false, py::arg("build_choice_origins") = false);
+    m.def("_build_sparse_model_from_drn", &storm::api::buildExplicitDRNModel<double>, "Build the model from DRN", py::arg("file"));
+    m.def("_build_sparse_parametric_model_from_drn", &storm::api::buildExplicitDRNModel<storm::RationalFunction>, "Build the parametric model from DRN", py::arg("file"));
+    m.def("build_sparse_model_from_explicit", &storm::api::buildExplicitModel<double>, "Build the model model from explicit input", py::arg("transition_file"), py::arg("labeling_file"), py::arg("state_reward_file") = "", py::arg("transition_reward_file") = "", py::arg("choice_labeling_file") = "");
 
 }
 
diff --git a/src/core/input.cpp b/src/core/input.cpp
index 67be172..b5152a3 100644
--- a/src/core/input.cpp
+++ b/src/core/input.cpp
@@ -15,7 +15,7 @@ void define_property(py::module& m) {
 void define_input(py::module& m) {
 
     // Parse prism program
-    m.def("parse_prism_program", &storm::parseProgram, "Parse prism program", py::arg("path"));
+    m.def("parse_prism_program", &storm::api::parseProgram, "Parse prism program", py::arg("path"));
 
     // PrismType
     py::enum_<storm::prism::Program::ModelType>(m, "PrismModelType", "Type of the prism model")
diff --git a/src/core/modelchecking.cpp b/src/core/modelchecking.cpp
index faa488c..f7d9d5e 100644
--- a/src/core/modelchecking.cpp
+++ b/src/core/modelchecking.cpp
@@ -1,14 +1,13 @@
 #include "modelchecking.h"
 #include "result.h"
 
-// Thin wrapper for model checking
-std::shared_ptr<storm::modelchecker::CheckResult> modelChecking(std::shared_ptr<storm::models::sparse::Model<double>> model, std::shared_ptr<storm::logic::Formula const> const& formula) {
-    return storm::verifySparseModel<double>(model, formula);
-}
+template<typename ValueType>
+using CheckTask = storm::modelchecker::CheckTask<storm::logic::Formula, ValueType>;
 
-// Thin wrapper for parametric model checking
-std::shared_ptr<storm::modelchecker::CheckResult> parametricModelChecking(std::shared_ptr<storm::models::sparse::Model<storm::RationalFunction>> model, std::shared_ptr<storm::logic::Formula const> const& formula) {
-    return storm::verifySparseModel<storm::RationalFunction>(model, formula);
+// Thin wrapper for model checking
+template<typename ValueType>
+std::shared_ptr<storm::modelchecker::CheckResult> modelCheckingSparseEngine(std::shared_ptr<storm::models::sparse::Model<ValueType>> model, CheckTask<ValueType> const& task) {
+    return storm::api::verifyWithSparseEngine<ValueType>(model, task);
 }
 
 // Thin wrapper for computing prob01 states
@@ -30,14 +29,23 @@ std::pair<storm::storage::BitVector, storm::storage::BitVector> computeProb01max
 // Define python bindings
 void define_modelchecking(py::module& m) {
 
+    // CheckTask
+    py::class_<CheckTask<double>, std::shared_ptr<CheckTask<double>>>(m, "CheckTask", "Task for model checking")
+    //m.def("create_check_task", &storm::api::createTask, "Create task for verification", py::arg("formula"), py::arg("only_initial_states") = false);
+        .def(py::init<storm::logic::Formula const&, bool>(), py::arg("formula"), py::arg("only_initial_states") = false)
+    ;
+    py::class_<CheckTask<storm::RationalFunction>, std::shared_ptr<CheckTask<storm::RationalFunction>>>(m, "ParametricCheckTask", "Task for parametric model checking")
+    //m.def("create_check_task", &storm::api::createTask, "Create task for verification", py::arg("formula"), py::arg("only_initial_states") = false);
+        .def(py::init<storm::logic::Formula const&, bool>(), py::arg("formula"), py::arg("only_initial_states") = false)
+    ;
+
     // Model checking
-    m.def("_model_checking", &modelChecking, "Perform model checking", py::arg("model"), py::arg("formula"));
-    m.def("_parametric_model_checking", &parametricModelChecking, "Perform parametric model checking", py::arg("model"), py::arg("formula"));
+    m.def("_model_checking_sparse_engine", &modelCheckingSparseEngine<double>, "Perform model checking", py::arg("model"), py::arg("task"));
+    m.def("_parametric_model_checking_sparse_engine", &modelCheckingSparseEngine<storm::RationalFunction>, "Perform parametric model checking", py::arg("model"), py::arg("task"));
     m.def("_compute_prob01states_double", &computeProb01<double>, "Compute prob-0-1 states", py::arg("model"), py::arg("phi_states"), py::arg("psi_states"));
     m.def("_compute_prob01states_rationalfunc", &computeProb01<storm::RationalFunction>, "Compute prob-0-1 states", py::arg("model"), py::arg("phi_states"), py::arg("psi_states"));
     m.def("_compute_prob01states_min_double", &computeProb01min<double>, "Compute prob-0-1 states (min)", py::arg("model"), py::arg("phi_states"), py::arg("psi_states"));
     m.def("_compute_prob01states_max_double", &computeProb01max<double>, "Compute prob-0-1 states (max)", py::arg("model"), py::arg("phi_states"), py::arg("psi_states"));
     m.def("_compute_prob01states_min_rationalfunc", &computeProb01min<storm::RationalFunction>, "Compute prob-0-1 states (min)", py::arg("model"), py::arg("phi_states"), py::arg("psi_states"));
     m.def("_compute_prob01states_max_rationalfunc", &computeProb01max<storm::RationalFunction>, "Compute prob-0-1 states (max)", py::arg("model"), py::arg("phi_states"), py::arg("psi_states"));
-    
 }
diff --git a/src/core/pla.cpp b/src/core/pla.cpp
index 5506f71..1f34f12 100644
--- a/src/core/pla.cpp
+++ b/src/core/pla.cpp
@@ -1,5 +1,6 @@
 #include "pla.h"
 #include "src/helpers.h"
+#include "storm/modelchecker/parametric/SparseDtmcRegionChecker.h"
 
 typedef storm::modelchecker::parametric::SparseDtmcRegionChecker<storm::models::sparse::Dtmc<storm::RationalFunction>, double, storm::RationalNumber> SparseDtmcRegionChecker;
 typedef storm::storage::ParameterRegion<storm::RationalFunction> Region;
@@ -25,8 +26,8 @@ std::set<storm::Polynomial> gatherDerivatives(storm::models::sparse::Dtmc<storm:
 
 // Define python bindings
 void define_pla(py::module& m) {
-    
-    // RegionCheckResult 
+
+    // RegionCheckResult
     py::enum_<storm::modelchecker::parametric::RegionCheckResult>(m, "RegionCheckResult", "Types of region check results")
         .value("EXISTSSAT", storm::modelchecker::parametric::RegionCheckResult::ExistsSat)
         .value("EXISTSVIOLATED", storm::modelchecker::parametric::RegionCheckResult::ExistsViolated)
diff --git a/tests/core/test_parse.py b/tests/core/test_parse.py
index b0332bf..0610c25 100644
--- a/tests/core/test_parse.py
+++ b/tests/core/test_parse.py
@@ -23,7 +23,7 @@ class TestParse:
         assert str(properties[0].raw_formula) == formula
     
     def test_parse_explicit_dtmc(self):
-        model = stormpy.parse_explicit_model(get_example_path("dtmc", "die.tra"), get_example_path("dtmc", "die.lab"))
+        model = stormpy.build_sparse_model_from_explicit(get_example_path("dtmc", "die.tra"), get_example_path("dtmc", "die.lab"))
         assert model.nr_states == 13
         assert model.nr_transitions == 20
         assert model.model_type == stormpy.ModelType.DTMC
@@ -31,7 +31,7 @@ class TestParse:
         assert type(model) is stormpy.SparseDtmc
     
     def test_parse_explicit_mdp(self):
-        model = stormpy.parse_explicit_model(get_example_path("mdp", "two_dice.tra"), get_example_path("mdp", "two_dice.lab"))
+        model = stormpy.build_sparse_model_from_explicit(get_example_path("mdp", "two_dice.tra"), get_example_path("mdp", "two_dice.lab"))
         assert model.nr_states == 169
         assert model.nr_transitions == 436
         assert model.model_type == stormpy.ModelType.MDP
diff --git a/tests/storage/test_matrix.py b/tests/storage/test_matrix.py
index c709214..07d222a 100644
--- a/tests/storage/test_matrix.py
+++ b/tests/storage/test_matrix.py
@@ -6,7 +6,7 @@ import math
 
 class TestMatrix:
     def test_sparse_matrix(self):
-        model = stormpy.parse_explicit_model(get_example_path("dtmc", "die.tra"), get_example_path("dtmc", "die.lab"))
+        model = stormpy.build_sparse_model_from_explicit(get_example_path("dtmc", "die.tra"), get_example_path("dtmc", "die.lab"))
         matrix = model.transition_matrix
         assert type(matrix) is stormpy.storage.SparseMatrix
         assert matrix.nr_rows == model.nr_states
@@ -17,7 +17,7 @@ class TestMatrix:
             assert e.value() == 0.5 or e.value() == 0 or (e.value() == 1 and e.column > 6)
 
     def test_backward_matrix(self):
-        model = stormpy.parse_explicit_model(get_example_path("dtmc", "die.tra"), get_example_path("dtmc", "die.lab"))
+        model = stormpy.build_sparse_model_from_explicit(get_example_path("dtmc", "die.tra"), get_example_path("dtmc", "die.lab"))
         matrix = model.backward_transition_matrix
         assert type(matrix) is stormpy.storage.SparseMatrix
         assert matrix.nr_rows == model.nr_states
@@ -28,7 +28,7 @@ class TestMatrix:
             assert e.value() == 0.5 or e.value() == 0 or (e.value() == 1 and e.column > 6)
 
     def test_change_sparse_matrix(self):
-        model = stormpy.parse_explicit_model(get_example_path("dtmc", "die.tra"), get_example_path("dtmc", "die.lab"))
+        model = stormpy.build_sparse_model_from_explicit(get_example_path("dtmc", "die.tra"), get_example_path("dtmc", "die.lab"))
         matrix = model.transition_matrix
         for e in matrix:
             assert e.value() == 0.5 or e.value() == 0 or e.value() == 1
@@ -43,7 +43,7 @@ class TestMatrix:
 
     def test_change_sparse_matrix_modelchecking(self):
         import stormpy.logic
-        model = stormpy.parse_explicit_model(get_example_path("dtmc", "die.tra"), get_example_path("dtmc", "die.lab"))
+        model = stormpy.build_sparse_model_from_explicit(get_example_path("dtmc", "die.tra"), get_example_path("dtmc", "die.lab"))
         matrix = model.transition_matrix
         # Check matrix
         for e in matrix:
diff --git a/tests/storage/test_state.py b/tests/storage/test_state.py
index 56375a8..9b64ef4 100644
--- a/tests/storage/test_state.py
+++ b/tests/storage/test_state.py
@@ -3,7 +3,7 @@ from helpers.helper import get_example_path
 
 class TestState:
     def test_states_dtmc(self):
-        model = stormpy.parse_explicit_model(get_example_path("dtmc", "die.tra"), get_example_path("dtmc", "die.lab"))
+        model = stormpy.build_sparse_model_from_explicit(get_example_path("dtmc", "die.tra"), get_example_path("dtmc", "die.lab"))
         i = 0
         states = model.states
         assert len(states) == 13
@@ -17,7 +17,7 @@ class TestState:
         assert state.id == 5
 
     def test_states_mdp(self):
-        model = stormpy.parse_explicit_model(get_example_path("mdp", "two_dice.tra"), get_example_path("mdp", "two_dice.lab"))
+        model = stormpy.build_sparse_model_from_explicit(get_example_path("mdp", "two_dice.tra"), get_example_path("mdp", "two_dice.lab"))
         i = 0
         assert len(model.states) == 169
         for state in model.states:
@@ -31,7 +31,7 @@ class TestState:
         assert state.id == 1
 
     def test_labels_dtmc(self):
-        model = stormpy.parse_explicit_model(get_example_path("dtmc", "die.tra"), get_example_path("dtmc", "die.lab"))
+        model = stormpy.build_sparse_model_from_explicit(get_example_path("dtmc", "die.tra"), get_example_path("dtmc", "die.lab"))
         labelsOrig = [["init"], [], [], [], [], [], [], ["one", "done"], ["two", "done"], ["three", "done"], ["four", "done"], ["five", "done"], ["six", "done"]]
         i = 0
         for state in model.states:
@@ -41,14 +41,14 @@ class TestState:
             i += 1
 
     def test_actions_dtmc(self):
-        model = stormpy.parse_explicit_model(get_example_path("dtmc", "die.tra"), get_example_path("dtmc", "die.lab"))
+        model = stormpy.build_sparse_model_from_explicit(get_example_path("dtmc", "die.tra"), get_example_path("dtmc", "die.lab"))
         for state in model.states:
             assert len(state.actions) == 1
             for action in state.actions:
                 assert action.id == 0
 
     def test_actions_mdp(self):
-        model = stormpy.parse_explicit_model(get_example_path("mdp", "two_dice.tra"), get_example_path("mdp", "two_dice.lab"))
+        model = stormpy.build_sparse_model_from_explicit(get_example_path("mdp", "two_dice.tra"), get_example_path("mdp", "two_dice.lab"))
         for state in model.states:
             assert len(state.actions) == 1 or len(state.actions) == 2
             for action in state.actions:
@@ -61,7 +61,7 @@ class TestState:
                             (6, 2, 0.5), (6, 12, 0.5), (7, 7, 1), (8, 8, 1),
                             (9, 9, 1), (10, 10, 1), (11, 11, 1), (12, 12, 1)
                            ]
-        model = stormpy.parse_explicit_model(get_example_path("dtmc", "die.tra"), get_example_path("dtmc", "die.lab"))
+        model = stormpy.build_sparse_model_from_explicit(get_example_path("dtmc", "die.tra"), get_example_path("dtmc", "die.lab"))
         i = 0
         for state in model.states:
             for action in state.actions:
@@ -74,7 +74,7 @@ class TestState:
                     assert transition.value() == transition_orig[2]
 
     def test_transitions_mdp(self):
-        model = stormpy.parse_explicit_model(get_example_path("mdp", "two_dice.tra"), get_example_path("mdp", "two_dice.lab"))
+        model = stormpy.build_sparse_model_from_explicit(get_example_path("mdp", "two_dice.tra"), get_example_path("mdp", "two_dice.lab"))
         assert model.states[1].id == 1
         for state in model.states:
             i = 0
@@ -91,7 +91,7 @@ class TestState:
                             (6, 2, 0.5), (6, 12, 0.5), (7, 7, 1), (8, 8, 1),
                             (9, 9, 1), (10, 10, 1), (11, 11, 1), (12, 12, 1)
                            ]
-        model = stormpy.parse_explicit_model(get_example_path("dtmc", "die.tra"), get_example_path("dtmc", "die.lab"))
+        model = stormpy.build_sparse_model_from_explicit(get_example_path("dtmc", "die.tra"), get_example_path("dtmc", "die.lab"))
         i = 0
         for state in model.states:
             for transition in model.transition_matrix.row_iter(state.id, state.id):