diff --git a/src/storm/transformer/ChoiceSelector.cpp b/src/storm/transformer/ChoiceSelector.cpp
index 1044cb42c..752cc1ee0 100644
--- a/src/storm/transformer/ChoiceSelector.cpp
+++ b/src/storm/transformer/ChoiceSelector.cpp
@@ -1,5 +1,8 @@
 #include "storm/transformer/ChoiceSelector.h"
 #include "storm/models/sparse/Mdp.h"
+#include "storm/models/sparse/MarkovAutomaton.h"
+
+#include "storm/exceptions/UnexpectedException.h"
 
 namespace  storm {
     namespace transformer {
@@ -18,7 +21,15 @@ namespace  storm {
             if (inputModel.hasChoiceOrigins()) {
                 newComponents.choiceOrigins = inputModel.getChoiceOrigins()->selectChoices(enabledActions);
             }
-            return std::make_shared<storm::models::sparse::Mdp<ValueType, RewardModelType>>(std::move(newComponents));
+            if (inputModel.isOfType(storm::models::ModelType::MarkovAutomaton)) {
+                auto const& ma = *inputModel.template as<storm::models::sparse::MarkovAutomaton<ValueType, RewardModelType>>();
+                newComponents.markovianStates = ma.getMarkovianStates();
+                newComponents.exitRates = ma.getExitRates();
+                return std::make_shared<storm::models::sparse::MarkovAutomaton<ValueType, RewardModelType>>(std::move(newComponents));
+            } else {
+                STORM_LOG_THROW(inputModel.isOfType(storm::models::ModelType::Mdp), storm::exceptions::UnexpectedException, "Unexpected model type for choice selector.");
+                return std::make_shared<storm::models::sparse::Mdp<ValueType, RewardModelType>>(std::move(newComponents));
+            }
         }
 
         template class ChoiceSelector<double>;
diff --git a/src/storm/transformer/SubsystemBuilder.cpp b/src/storm/transformer/SubsystemBuilder.cpp
new file mode 100644
index 000000000..b8ba5f0e3
--- /dev/null
+++ b/src/storm/transformer/SubsystemBuilder.cpp
@@ -0,0 +1,127 @@
+#include "storm/transformer/SubsystemBuilder.h"
+
+#include <boost/optional.hpp>
+#include <storm/exceptions/UnexpectedException.h>
+
+#include "storm/models/sparse/StandardRewardModel.h"
+#include "storm/utility/constants.h"
+#include "storm/utility/graph.h"
+#include "storm/utility/macros.h"
+#include "storm/utility/vector.h"
+#include "storm/models/sparse/Dtmc.h"
+#include "storm/models/sparse/Mdp.h"
+#include "storm/models/sparse/Ctmc.h"
+#include "storm/models/sparse/MarkovAutomaton.h"
+#include "storm/storage/sparse/ModelComponents.h"
+#include "storm/utility/builder.h"
+
+#include "storm/exceptions/InvalidArgumentException.h"
+#include "storm/exceptions/UnexpectedException.h"
+
+namespace storm {
+    namespace transformer {
+        
+        template <typename ValueType, typename RewardModelType>
+        void transformModelSpecificComponents(storm::models::sparse::Model<ValueType, RewardModelType> const& originalModel,
+                               storm::storage::BitVector const& subsystem,
+                               storm::storage::sparse::ModelComponents<ValueType, RewardModelType>& components) {
+            if (originalModel.isOfType(storm::models::ModelType::MarkovAutomaton)) {
+                auto const& ma = *originalModel.template as<storm::models::sparse::MarkovAutomaton<ValueType, RewardModelType>>();
+                components.markovianStates = ma.getMarkovianStates() % subsystem;
+                components.exitRates = storm::utility::vector::filterVector(ma.getExitRates(), subsystem);
+                components.rateTransitions = false; // Note that originalModel.getTransitionMatrix() contains probabilities
+            } else if (originalModel.isOfType(storm::models::ModelType::Ctmc)) {
+                auto const& ctmc = *originalModel.template as<storm::models::sparse::Ctmc<ValueType, RewardModelType>>();
+                components.exitRates = storm::utility::vector::filterVector(ctmc.getExitRateVector(), subsystem);
+                components.rateTransitions = true;
+            } else {
+                STORM_LOG_THROW(originalModel.isOfType(storm::models::ModelType::Dtmc) || originalModel.isOfType(storm::models::ModelType::Mdp), storm::exceptions::UnexpectedException, "Unexpected model type.");
+            }
+        }
+        
+        template<typename RewardModelType>
+        RewardModelType transformRewardModel(RewardModelType const& originalRewardModel, storm::storage::BitVector const& subsystem, storm::storage::BitVector const& subsystemActions) {
+            boost::optional<std::vector<typename RewardModelType::ValueType>> stateRewardVector;
+            boost::optional<std::vector<typename RewardModelType::ValueType>> stateActionRewardVector;
+            boost::optional<storm::storage::SparseMatrix<typename RewardModelType::ValueType>> transitionRewardMatrix;
+            if (originalRewardModel.hasStateRewards()){
+                stateRewardVector = storm::utility::vector::filterVector(originalRewardModel.getStateRewardVector(), subsystem);
+            }
+            if (originalRewardModel.hasStateActionRewards()){
+                stateActionRewardVector = storm::utility::vector::filterVector(originalRewardModel.getStateActionRewardVector(), subsystemActions);
+            }
+            if (originalRewardModel.hasTransitionRewards()){
+                transitionRewardMatrix = originalRewardModel.getTransitionRewardMatrix().getSubmatrix(false, subsystemActions, subsystem);
+            }
+            return RewardModelType(std::move(stateRewardVector), std::move(stateActionRewardVector), std::move(transitionRewardMatrix));
+        }
+
+        template <typename ValueType, typename RewardModelType>
+        SubsystemBuilderReturnType<ValueType, RewardModelType> internalBuildSubsystem(storm::models::sparse::Model<ValueType, RewardModelType> const& originalModel, storm::storage::BitVector const& subsystemStates, storm::storage::BitVector const& subsystemActions) {
+
+            SubsystemBuilderReturnType<ValueType, RewardModelType> result;
+            uint_fast64_t subsystemStateCount = subsystemStates.getNumberOfSetBits();
+            result.newToOldStateIndexMapping.reserve(subsystemStateCount);
+            result.keptActions = storm::storage::BitVector(originalModel.getTransitionMatrix().getRowCount(), false);
+            for (auto subsysState : subsystemStates) {
+                result.newToOldStateIndexMapping.push_back(subsysState);
+                bool stateHasOneChoiceLeft = false;
+                for (uint_fast64_t row = subsystemActions.getNextSetIndex(originalModel.getTransitionMatrix().getRowGroupIndices()[subsysState]); row < originalModel.getTransitionMatrix().getRowGroupIndices()[subsysState+1]; row = subsystemActions.getNextSetIndex(row+1)) {
+                    bool allRowEntriesStayInSubsys = true;
+                    for (auto const& entry : originalModel.getTransitionMatrix().getRow(row)) {
+                        if (!subsystemStates.get(entry.getColumn())) {
+                            allRowEntriesStayInSubsys = false;
+                            break;
+                        }
+                    }
+                    stateHasOneChoiceLeft |= allRowEntriesStayInSubsys;
+                    result.keptActions.set(row, allRowEntriesStayInSubsys);
+                }
+                 STORM_LOG_THROW(stateHasOneChoiceLeft, storm::exceptions::InvalidArgumentException, "The subsystem would contain a deadlock state.");
+            }
+            
+            // Transform the components of the model
+            storm::storage::sparse::ModelComponents<ValueType, RewardModelType> components;
+            components.transitionMatrix = originalModel.getTransitionMatrix().getSubmatrix(false, result.keptActions, subsystemStates);
+            components.stateLabeling = originalModel.getStateLabeling().getSubLabeling(subsystemStates);
+            for (auto const& rewardModel : originalModel.getRewardModels()){
+                components.rewardModels.insert(std::make_pair(rewardModel.first, transformRewardModel(rewardModel.second, subsystemStates, result.keptActions)));
+            }
+            if (originalModel.hasChoiceLabeling()) {
+                components.choiceLabeling = originalModel.getChoiceLabeling().getSubLabeling(result.keptActions);
+            }
+            if (originalModel.hasStateValuations()) {
+                components.stateValuations = originalModel.getStateValuations().selectStates(subsystemStates);
+            }
+            if (originalModel.hasChoiceOrigins()) {
+                components.choiceOrigins = originalModel.getChoiceOrigins()->selectChoices(result.keptActions);
+            }
+            
+            transformModelSpecificComponents<ValueType, RewardModelType>(originalModel, subsystemStates, components);
+            
+            result.model = storm::utility::builder::buildModelFromComponents(originalModel.getType(), std::move(components));
+            STORM_LOG_DEBUG("Subsystem Builder is done. Resulting model has " << result.model->getNumberOfStates() << " states.");
+            return result;
+        }
+        
+        template <typename ValueType, typename RewardModelType>
+        SubsystemBuilderReturnType<ValueType, RewardModelType> buildSubsystem(storm::models::sparse::Model<ValueType, RewardModelType> const& originalModel, storm::storage::BitVector const& subsystemStates, storm::storage::BitVector const& subsystemActions, bool keepUnreachableStates) {
+            STORM_LOG_DEBUG("Invoked subsystem builder on model with " << originalModel.getNumberOfStates() << " states.");
+            storm::storage::BitVector initialStates = originalModel.getInitialStates() & subsystemStates;
+            STORM_LOG_THROW(!initialStates.empty(), storm::exceptions::InvalidArgumentException, "The subsystem would not contain any initial states");
+            
+            STORM_LOG_THROW(!subsystemStates.empty(), storm::exceptions::InvalidArgumentException, "Invoked SubsystemBuilder for an empty subsystem.");
+            if (keepUnreachableStates) {
+                return internalBuildSubsystem(originalModel, subsystemStates, subsystemActions);
+            } else {
+                auto actualSubsystem = storm::utility::graph::getReachableStates(originalModel.getTransitionMatrix(), initialStates, subsystemStates, storm::storage::BitVector(subsystemStates.size(), false), false, 0, subsystemActions);
+                return internalBuildSubsystem(originalModel, actualSubsystem, subsystemActions);
+            }
+        }
+        
+        template SubsystemBuilderReturnType<double> buildSubsystem(storm::models::sparse::Model<double> const& originalModel, storm::storage::BitVector const& subsystemStates, storm::storage::BitVector const& subsystemActions, bool keepUnreachableStates = true);
+        template SubsystemBuilderReturnType<double, storm::models::sparse::StandardRewardModel<storm::Interval>> buildSubsystem(storm::models::sparse::Model<double, storm::models::sparse::StandardRewardModel<storm::Interval>> const& originalModel, storm::storage::BitVector const& subsystemStates, storm::storage::BitVector const& subsystemActions, bool keepUnreachableStates = true);
+        template SubsystemBuilderReturnType<storm::RationalNumber> buildSubsystem(storm::models::sparse::Model<storm::RationalNumber> const& originalModel, storm::storage::BitVector const& subsystemStates, storm::storage::BitVector const& subsystemActions, bool keepUnreachableStates = true);
+        template SubsystemBuilderReturnType<storm::RationalFunction> buildSubsystem(storm::models::sparse::Model<storm::RationalFunction> const& originalModel, storm::storage::BitVector const& subsystemStates, storm::storage::BitVector const& subsystemActions, bool keepUnreachableStates = true);
+    }
+}
diff --git a/src/storm/transformer/SubsystemBuilder.h b/src/storm/transformer/SubsystemBuilder.h
index 1b5e7bd53..cc6431745 100644
--- a/src/storm/transformer/SubsystemBuilder.h
+++ b/src/storm/transformer/SubsystemBuilder.h
@@ -1,170 +1,43 @@
-#ifndef STORM_TRANSFORMER_SUBSYSTEMBUILDER_H
-#define STORM_TRANSFORMER_SUBSYSTEMBUILDER_H
-
+#pragma once
 
 #include <memory>
-#include <boost/optional.hpp>
+#include <vector>
 
+#include "storm/storage/BitVector.h"
+#include "storm/models/sparse/Model.h"
 #include "storm/models/sparse/StandardRewardModel.h"
-#include "storm/utility/constants.h"
-#include "storm/utility/graph.h"
-#include "storm/utility/macros.h"
-#include "storm/utility/vector.h"
-#include "storm/models/sparse/Dtmc.h"
-#include "storm/models/sparse/Mdp.h"
-#include "storm/models/sparse/Ctmc.h"
-#include "storm/models/sparse/MarkovAutomaton.h"
-#include "storm/storage/sparse/ModelComponents.h"
-#include "storm/utility/builder.h"
-
-#include "storm/exceptions/InvalidArgumentException.h"
-#include "storm/exceptions/InvalidStateException.h"
 
 namespace storm {
     namespace transformer {
         
+        template <typename ValueType, typename RewardModelType = storm::models::sparse::StandardRewardModel<ValueType>>
+        struct SubsystemBuilderReturnType {
+            // The resulting model
+            std::shared_ptr<storm::models::sparse::Model<ValueType, RewardModelType>> model;
+            // Gives for each state in the resulting model the corresponding state in the original model.
+            std::vector<uint_fast64_t> newToOldStateIndexMapping;
+            // marks the actions of the original model that are still available in the subsystem
+            storm::storage::BitVector keptActions;
+        };
+        
         /*
-         * Removes all states that are not part of the subsystem
+         * Removes all states and actions that are not part of the subsystem.
+         * A state is part of the subsystem iff
+         *    * it is selected in subsystemStates AND
+         *    * keepUnreachableStates is true or the state is reachable from the initial states
+         * An action is part of the subsystem iff
+         *    * it is selected in subsystemActions AND
+         *    * it originates from a state that is part of the subsystem AND
+         *    * it does not contain a transition leading to a state outside of the subsystem.
+         *
+         * If this introduces a deadlock state (i.e., a state without an action) an exception is thrown.
+         *
+         * @param originalModel The original model.
+         * @param subsystemStates The selected states.
+         * @param subsystemActions The selected actions
+         * @param keepUnreachableStates if true, states that are not reachable from the initial state are kept
          */
-        template <typename SparseModelType>
-        class SubsystemBuilder {
-        public:
-
-            struct SubsystemBuilderReturnType {
-                // The resulting model
-                std::shared_ptr<SparseModelType> model;
-                // Gives for each state in the resulting model the corresponding state in the original model.
-                std::vector<uint_fast64_t> newToOldStateIndexMapping;
-                // marks the actions of the original model that are still available in the subsystem
-                storm::storage::BitVector keptActions;
-            };
-            
-            /*
-             * Removes all states and actions that are not part of the subsystem.
-             * A state is part of the subsystem iff it is selected in subsystemStates.
-             * An action is part of the subsystem iff 
-             *    * it is selected in subsystemActions AND
-             *    * it originates from a state that is part of the subsystem AND
-             *    * it does not contain a transition leading to a state outside of the subsystem.
-             *
-             * If this introduces a deadlock state (i.e., a state without an action) an exception is thrown.
-             * 
-             * @param originalModel The original model.
-             * @param subsystemStates The selected states.
-             * @param subsystemActions The selected actions
-             */
-            static SubsystemBuilderReturnType transform(SparseModelType const& originalModel, storm::storage::BitVector const& subsystemStates, storm::storage::BitVector const& subsystemActions) {
-                STORM_LOG_DEBUG("Invoked subsystem builder on model with " << originalModel.getNumberOfStates() << " states.");
-                STORM_LOG_THROW(!(originalModel.getInitialStates() & subsystemStates).empty(), storm::exceptions::InvalidArgumentException, "The subsystem would not contain any initial states");
-                SubsystemBuilderReturnType result;
-                
-                uint_fast64_t subsystemStateCount = subsystemStates.getNumberOfSetBits();
-                STORM_LOG_THROW(subsystemStateCount != 0, storm::exceptions::InvalidArgumentException, "Invoked SubsystemBuilder for an empty subsystem.");
-                if (subsystemStateCount == subsystemStates.size() && subsystemActions.full()) {
-                    result.model = std::make_shared<SparseModelType>(originalModel);
-                    result.newToOldStateIndexMapping = storm::utility::vector::buildVectorForRange(0, result.model->getNumberOfStates());
-                    result.keptActions = storm::storage::BitVector(result.model->getTransitionMatrix().getRowCount(), true);
-                    return result;
-                }
-                
-                result.newToOldStateIndexMapping.reserve(subsystemStateCount);
-                result.keptActions = storm::storage::BitVector(originalModel.getTransitionMatrix().getRowCount(), false);
-                for (auto subsysState : subsystemStates) {
-                    result.newToOldStateIndexMapping.push_back(subsysState);
-                    bool stateHasOneChoiceLeft = false;
-                    for (uint_fast64_t row = subsystemActions.getNextSetIndex(originalModel.getTransitionMatrix().getRowGroupIndices()[subsysState]); row < originalModel.getTransitionMatrix().getRowGroupIndices()[subsysState+1]; row = subsystemActions.getNextSetIndex(row+1)) {
-                        bool allRowEntriesStayInSubsys = true;
-                        for (auto const& entry : originalModel.getTransitionMatrix().getRow(row)) {
-                            if (!subsystemStates.get(entry.getColumn())) {
-                                allRowEntriesStayInSubsys = false;
-                                break;
-                            }
-                        }
-                        stateHasOneChoiceLeft |= allRowEntriesStayInSubsys;
-                        result.keptActions.set(row, allRowEntriesStayInSubsys);
-                    }
-                     STORM_LOG_THROW(stateHasOneChoiceLeft, storm::exceptions::InvalidArgumentException, "The subsystem would contain a deadlock state.");
-                }
-                
-                // Transform the components of the model
-                storm::storage::sparse::ModelComponents<typename SparseModelType::ValueType, typename SparseModelType::RewardModelType> components;
-                components.transitionMatrix = originalModel.getTransitionMatrix().getSubmatrix(false, result.keptActions, subsystemStates);
-                components.stateLabeling = originalModel.getStateLabeling().getSubLabeling(subsystemStates);
-                for (auto const& rewardModel : originalModel.getRewardModels()){
-                    components.rewardModels.insert(std::make_pair(rewardModel.first, transformRewardModel(rewardModel.second, subsystemStates, result.keptActions)));
-                }
-                if (originalModel.hasChoiceLabeling()) {
-                    components.choiceLabeling = originalModel.getChoiceLabeling().getSubLabeling(result.keptActions);
-                }
-                if (originalModel.hasStateValuations()) {
-                    components.stateValuations = originalModel.getStateValuations().selectStates(subsystemStateCount);
-                }
-                if (originalModel.hasChoiceOrigins()) {
-                    components.choiceOrigins = originalModel.getChoiceOrigins().selectChoices(result.keptActions);
-                }
-                
-                transformModelSpecificComponents(originalModel, subsystemStates, components);
-                
-                result.model = storm::utility::builder::buildModelFromComponents(originalModel.getType(), std::move(components));
-                STORM_LOG_DEBUG("Subsystem Builder is done. Resulting model has " << result.model->getNumberOfStates() << " states.");
-                return result;
-            }
-            
-        private:
-            template<typename ValueType = typename SparseModelType::RewardModelType::ValueType, typename RewardModelType = typename SparseModelType::RewardModelType>
-            static RewardModelType transformRewardModel(RewardModelType const& originalRewardModel, storm::storage::BitVector const& subsystem, storm::storage::BitVector const& subsystemActions) {
-                boost::optional<std::vector<ValueType>> stateRewardVector;
-                boost::optional<std::vector<ValueType>> stateActionRewardVector;
-                boost::optional<storm::storage::SparseMatrix<ValueType>> transitionRewardMatrix;
-                if (originalRewardModel.hasStateRewards()){
-                    stateRewardVector = storm::utility::vector::filterVector(originalRewardModel.getStateRewardVector(), subsystem);
-                }
-                if (originalRewardModel.hasStateActionRewards()){
-                    stateActionRewardVector = storm::utility::vector::filterVector(originalRewardModel.getStateActionRewardVector(), subsystemActions);
-                }
-                if (originalRewardModel.hasTransitionRewards()){
-                    transitionRewardMatrix = originalRewardModel.getTransitionRewardMatrix().getSubmatrix(false, subsystemActions, subsystem);
-                }
-                return RewardModelType(std::move(stateRewardVector), std::move(stateActionRewardVector), std::move(transitionRewardMatrix));
-            }
-            
-            template<typename MT = SparseModelType>
-            static typename std::enable_if<
-            std::is_same<MT,storm::models::sparse::Dtmc<typename SparseModelType::ValueType>>::value ||
-            std::is_same<MT,storm::models::sparse::Mdp<typename SparseModelType::ValueType>>::value,
-            MT>::type
-            transformModelSpecificComponents(MT const&,
-                                   storm::storage::BitVector const& /*subsystem*/,
-                                   storm::storage::sparse::ModelComponents<typename SparseModelType::ValueType, typename SparseModelType::RewardModelType>&) {
-                // Intentionally left empty
-            }
-            
-            template<typename MT = SparseModelType>
-            static typename std::enable_if<
-            std::is_same<MT,storm::models::sparse::Ctmc<typename SparseModelType::ValueType>>::value,
-            MT>::type
-            transformModelSpecificComponents(MT const& originalModel,
-                                   storm::storage::BitVector const& subsystem,
-                                   storm::storage::SparseMatrix<typename MT::ValueType>& matrix,
-                                   storm::storage::sparse::ModelComponents<typename SparseModelType::ValueType, typename SparseModelType::RewardModelType>& components) {
-                components.exitRates = storm::utility::vector::filterVector(originalModel.getExitRateVector(), subsystem);
-                components.rateTransitions = true;
-            }
-        
-            template<typename MT = SparseModelType>
-            static typename std::enable_if<
-            std::is_same<MT,storm::models::sparse::MarkovAutomaton<typename SparseModelType::ValueType>>::value,
-            MT>::type
-            transformModelSpecificComponents(MT const& originalModel,
-                                   storm::storage::BitVector const& subsystem,
-                                   storm::storage::SparseMatrix<typename MT::ValueType>& matrix,
-                                   storm::storage::sparse::ModelComponents<typename SparseModelType::ValueType, typename SparseModelType::RewardModelType>& components) {
-                components.markovianStates = originalModel.getMarkovianStates() % subsystem;
-                components.exitRates = storm::utility::vector::filterVector(originalModel.getExitRates(), subsystem);
-                components.rateTransitions = false; // Note that originalModel.getTransitionMatrix() contains probabilities
-            }
-            
-        };
+        template <typename ValueType, typename RewardModelType = storm::models::sparse::StandardRewardModel<ValueType>>
+        SubsystemBuilderReturnType<ValueType, RewardModelType> buildSubsystem(storm::models::sparse::Model<ValueType, RewardModelType> const& originalModel, storm::storage::BitVector const& subsystemStates, storm::storage::BitVector const& subsystemActions, bool keepUnreachableStates = true);
     }
 }
-#endif // STORM_TRANSFORMER_SUBSYSTEMBUILDER_H
diff --git a/src/storm/utility/graph.cpp b/src/storm/utility/graph.cpp
index 5bb8cc758..076185012 100644
--- a/src/storm/utility/graph.cpp
+++ b/src/storm/utility/graph.cpp
@@ -30,7 +30,7 @@ namespace storm {
         namespace graph {
             
             template<typename T>
-            storm::storage::BitVector getReachableStates(storm::storage::SparseMatrix<T> const& transitionMatrix, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& constraintStates, storm::storage::BitVector const& targetStates, bool useStepBound, uint_fast64_t maximalSteps) {
+            storm::storage::BitVector getReachableStates(storm::storage::SparseMatrix<T> const& transitionMatrix, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& constraintStates, storm::storage::BitVector const& targetStates, bool useStepBound, uint_fast64_t maximalSteps, boost::optional<storm::storage::BitVector> const& choiceFilter) {
                 storm::storage::BitVector reachableStates(initialStates);
                 
                 uint_fast64_t numberOfStates = transitionMatrix.getRowGroupCount();
@@ -64,26 +64,38 @@ namespace storm {
                             continue;
                         }
                     }
-                
-                    for (auto const& successor : transitionMatrix.getRowGroup(currentState)) {
-                        // Only explore the state if the transition was actually there and the successor has not yet
-                        // been visited.
-                        if (!storm::utility::isZero(successor.getValue()) && (!reachableStates.get(successor.getColumn()) || (useStepBound && remainingSteps[successor.getColumn()] < currentStepBound - 1))) {
-                            // If the successor is one of the target states, we need to include it, but must not explore
-                            // it further.
-                            if (targetStates.get(successor.getColumn())) {
-                                reachableStates.set(successor.getColumn());
-                            } else if (constraintStates.get(successor.getColumn())) {
-                                // However, if the state is in the constrained set of states, we potentially need to follow it.
-                                if (useStepBound) {
-                                    // As there is at least one more step to go, we need to push the state and the new number of steps.
-                                    remainingSteps[successor.getColumn()] = currentStepBound - 1;
-                                    stepStack.push_back(currentStepBound - 1);
+                    
+                    
+                    uint64_t row = transitionMatrix.getRowGroupIndices()[currentState];
+                    if (choiceFilter) {
+                        row = choiceFilter->getNextSetIndex(row);
+                    }
+                    uint64_t const rowGroupEnd = transitionMatrix.getRowGroupIndices()[currentState + 1];
+                    while (row < rowGroupEnd) {
+                        for (auto const& successor : transitionMatrix.getRow(row)) {
+                            // Only explore the state if the transition was actually there and the successor has not yet
+                            // been visited.
+                            if (!storm::utility::isZero(successor.getValue()) && (!reachableStates.get(successor.getColumn()) || (useStepBound && remainingSteps[successor.getColumn()] < currentStepBound - 1))) {
+                                // If the successor is one of the target states, we need to include it, but must not explore
+                                // it further.
+                                if (targetStates.get(successor.getColumn())) {
+                                    reachableStates.set(successor.getColumn());
+                                } else if (constraintStates.get(successor.getColumn())) {
+                                    // However, if the state is in the constrained set of states, we potentially need to follow it.
+                                    if (useStepBound) {
+                                        // As there is at least one more step to go, we need to push the state and the new number of steps.
+                                        remainingSteps[successor.getColumn()] = currentStepBound - 1;
+                                        stepStack.push_back(currentStepBound - 1);
+                                    }
+                                    reachableStates.set(successor.getColumn());
+                                    stack.push_back(successor.getColumn());
                                 }
-                                reachableStates.set(successor.getColumn());
-                                stack.push_back(successor.getColumn());
                             }
                         }
+                        ++row;
+                        if (choiceFilter) {
+                            row = choiceFilter->getNextSetIndex(row);
+                        }
                     }
                 }
                 
@@ -1357,7 +1369,7 @@ namespace storm {
             }
 
 
-            template storm::storage::BitVector getReachableStates(storm::storage::SparseMatrix<double> const& transitionMatrix, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& constraintStates, storm::storage::BitVector const& targetStates, bool useStepBound, uint_fast64_t maximalSteps);
+            template storm::storage::BitVector getReachableStates(storm::storage::SparseMatrix<double> const& transitionMatrix, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& constraintStates, storm::storage::BitVector const& targetStates, bool useStepBound, uint_fast64_t maximalSteps, boost::optional<storm::storage::BitVector> const& choiceFilter);
             
             template storm::storage::BitVector getBsccCover(storm::storage::SparseMatrix<double> const& transitionMatrix);
            
@@ -1428,7 +1440,7 @@ namespace storm {
             
             // Instantiations for storm::RationalNumber.
 #ifdef STORM_HAVE_CARL
-            template storm::storage::BitVector getReachableStates(storm::storage::SparseMatrix<storm::RationalNumber> const& transitionMatrix, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& constraintStates, storm::storage::BitVector const& targetStates, bool useStepBound, uint_fast64_t maximalSteps);
+            template storm::storage::BitVector getReachableStates(storm::storage::SparseMatrix<storm::RationalNumber> const& transitionMatrix, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& constraintStates, storm::storage::BitVector const& targetStates, bool useStepBound, uint_fast64_t maximalSteps, boost::optional<storm::storage::BitVector> const& choiceFilter);
             
             template storm::storage::BitVector getBsccCover(storm::storage::SparseMatrix<storm::RationalNumber> const& transitionMatrix);
            
@@ -1479,7 +1491,7 @@ namespace storm {
             template std::vector<uint_fast64_t> getTopologicalSort(storm::storage::SparseMatrix<storm::RationalNumber> const& matrix);
             // End of instantiations for storm::RationalNumber.
             
-            template storm::storage::BitVector getReachableStates(storm::storage::SparseMatrix<storm::RationalFunction> const& transitionMatrix, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& constraintStates, storm::storage::BitVector const& targetStates, bool useStepBound, uint_fast64_t maximalSteps);
+            template storm::storage::BitVector getReachableStates(storm::storage::SparseMatrix<storm::RationalFunction> const& transitionMatrix, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& constraintStates, storm::storage::BitVector const& targetStates, bool useStepBound, uint_fast64_t maximalSteps, boost::optional<storm::storage::BitVector> const& choiceFilter);
             
             template storm::storage::BitVector getBsccCover(storm::storage::SparseMatrix<storm::RationalFunction> const& transitionMatrix);
             
diff --git a/src/storm/utility/graph.h b/src/storm/utility/graph.h
index 79103700b..1f7d4394b 100644
--- a/src/storm/utility/graph.h
+++ b/src/storm/utility/graph.h
@@ -62,9 +62,10 @@ namespace storm {
              * @param targetStates The target states that may not be passed.
              * @param useStepBound A flag that indicates whether or not to use the given number of maximal steps for the search.
              * @param maximalSteps The maximal number of steps to reach the psi states.
+             * @param choiceFilter If given, only choices for which the bitvector is true are considered.
              */
             template<typename T>
-            storm::storage::BitVector getReachableStates(storm::storage::SparseMatrix<T> const& transitionMatrix, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& constraintStates, storm::storage::BitVector const& targetStates, bool useStepBound = false, uint_fast64_t maximalSteps = 0);
+            storm::storage::BitVector getReachableStates(storm::storage::SparseMatrix<T> const& transitionMatrix, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& constraintStates, storm::storage::BitVector const& targetStates, bool useStepBound = false, uint_fast64_t maximalSteps = 0, boost::optional<storm::storage::BitVector> const& choiceFilter = boost::none);
 
             /*!
              * Retrieves a set of states that covers als BSCCs of the system in the sense that for every BSCC exactly