diff --git a/source/rrRoadRunner.cpp b/source/rrRoadRunner.cpp
index e26ec034ef..33104c041d 100644
--- a/source/rrRoadRunner.cpp
+++ b/source/rrRoadRunner.cpp
@@ -229,6 +229,11 @@ namespace rr {
 
         std::vector<SelectionRecord> mSelectionList;
 
+        /**
+         * the reaction rates of an output row, see getSelectedReactionRates
+         */
+        std::vector<double> selectedReactionRates;
+
         /**
          * ModelGenerator obtained from the factory
          */
@@ -1410,16 +1415,31 @@ namespace rr {
         return dResult;
     }
 
-    double RoadRunner::getNthSelectedOutput(size_t index, double currentTime) {
+    double RoadRunner::getNthSelectedOutput(size_t index, double currentTime,
+                                            const double *reactionRates) {
         const SelectionRecord &record = impl->mSelectionList[index];
 
         if (record.selectionType == SelectionRecord::TIME) {
             return currentTime;
+        } else if (reactionRates && record.selectionType == SelectionRecord::REACTION_RATE) {
+            return reactionRates[record.index];
         } else {
             return getValue(record);
         }
     }
 
+    const double *RoadRunner::getSelectedReactionRates() {
+        for (const SelectionRecord &record: impl->mSelectionList) {
+            if (record.selectionType == SelectionRecord::REACTION_RATE) {
+                std::vector<double> &rates = impl->selectedReactionRates;
+                rates.resize(impl->model->getNumReactions());
+                impl->model->getReactionRates(rates.size(), nullptr, rates.data());
+                return rates.data();
+            }
+        }
+        return nullptr;
+    }
+
     int RoadRunner::getTimeRowIndex() {
         for (u_int j = 0; j < impl->mSelectionList.size(); j++) {
             const SelectionRecord &record = impl->mSelectionList[j];
@@ -1432,8 +1452,9 @@ namespace rr {
     }
 
     void RoadRunner::getSelectedValues(ls::DoubleMatrix &results, int nRow, double currentTime) {
+        const double *reactionRates = getSelectedReactionRates();
         for (u_int j = 0; j < impl->mSelectionList.size(); j++) {
-            double out = getNthSelectedOutput(j, currentTime);
+            double out = getNthSelectedOutput(j, currentTime, reactionRates);
             results(nRow, j) = out;
         }
     }
@@ -1443,9 +1464,10 @@ namespace rr {
         assert(results.size() == impl->mSelectionList.size()
                && "given std::vector and selection list different size");
 
+        const double *reactionRates = getSelectedReactionRates();
         size_t size = results.size();
         for (size_t i = 0; i < size; ++i) {
-            results[i] = getNthSelectedOutput(i, currentTime);
+            results[i] = getNthSelectedOutput(i, currentTime, reactionRates);
         }
     }
 
@@ -5222,8 +5244,9 @@ namespace rr {
         std::vector<double> result;
         result.resize(impl->mSelectionList.size());
 
+        const double *reactionRates = getSelectedReactionRates();
         for (int i = 0; i < impl->mSelectionList.size(); i++) {
-            result[i] = getNthSelectedOutput(i, impl->model->getTime());
+            result[i] = getNthSelectedOutput(i, impl->model->getTime(), reactionRates);
         }
         return result;
     }
diff --git a/source/rrRoadRunner.h b/source/rrRoadRunner.h
index 3fb94a6ab4..a847aa00e9 100644
--- a/source/rrRoadRunner.h
+++ b/source/rrRoadRunner.h
@@ -1906,7 +1906,21 @@ namespace rr {
         bool populateResult();
 
 
-        double getNthSelectedOutput(size_t index, double currentTime);
+        /**
+         * the value of the n'th selection, a reaction rate from reactionRates
+         * if given, see getSelectedReactionRates
+         */
+        double getNthSelectedOutput(size_t index, double currentTime,
+                                    const double *reactionRates = nullptr);
+
+        /**
+         * Evaluates the reaction rates once for the selections of an output
+         * row and returns them, or nullptr if no reaction rate is selected.
+         * The rate of one reaction is evaluated with all reaction rates of the
+         * model, so reading every selected rate on its own made an output row
+         * quadratic in the number of reactions.
+         */
+        const double *getSelectedReactionRates();
 
         bool isParameterUsed(const std::string &sid);
 
diff --git a/test/cxx_api_tests/RoadRunnerAPITests.h b/test/cxx_api_tests/RoadRunnerAPITests.h
index d79c8801b8..a6b91a88a9 100644
--- a/test/cxx_api_tests/RoadRunnerAPITests.h
+++ b/test/cxx_api_tests/RoadRunnerAPITests.h
@@ -588,6 +588,31 @@ class RoadRunnerAPITests : public RoadRunnerTest {
         delete testModel;
     }
 
+    void getSelectedReactionRates() {
+        // the reaction rates of every output row are the ones of its state,
+        // also selected twice and between other selections
+        TestModel *testModel = TestModelFactory("SimpleFlux");
+        RoadRunner rr(testModel->str());
+        rr.setSelections({"time", "_J1", "[S1]", "_J0", "[S2]", "_J1"});
+        const ls::DoubleMatrix *result = rr.simulate(0, 10, 11);
+        ASSERT_EQ(result->numRows(), 11);
+        for (int i = 0; i < result->numRows(); ++i) {
+            const double s1 = (*result)(i, 2);
+            const double s2 = (*result)(i, 4);
+            EXPECT_DOUBLE_EQ((*result)(i, 1), 0.01 * s2);
+            EXPECT_DOUBLE_EQ((*result)(i, 3), 0.1 * s1);
+            EXPECT_DOUBLE_EQ((*result)(i, 5), 0.01 * s2);
+        }
+        // the last row is the state of the model at the end
+        std::vector<double> x = rr.getSelectedValues();
+        ASSERT_EQ(x.size(), 6);
+        EXPECT_EQ(x[1], rr.getValue("_J1"));
+        EXPECT_EQ(x[3], rr.getValue("_J0"));
+        EXPECT_EQ(x[5], rr.getValue("_J1"));
+        EXPECT_EQ(x[3], (*result)(10, 3));
+        delete testModel;
+    }
+
     void getIndependentFloatingSpeciesIds() {
         TestModel *testModel = TestModelFactory("SimpleFlux");
         RoadRunner rr(testModel->str());
diff --git a/test/cxx_api_tests/RoadRunnerAPITestsWithLLJit.cpp b/test/cxx_api_tests/RoadRunnerAPITestsWithLLJit.cpp
index 343f2b2659..f06e109093 100644
--- a/test/cxx_api_tests/RoadRunnerAPITestsWithLLJit.cpp
+++ b/test/cxx_api_tests/RoadRunnerAPITestsWithLLJit.cpp
@@ -194,6 +194,10 @@ TEST_F(RoadRunnerAPITestsWithLLJit, getSelectedValues){
     getSelectedValues();
 }
 
+TEST_F(RoadRunnerAPITestsWithLLJit, getSelectedReactionRates){
+    getSelectedReactionRates();
+}
+
 TEST_F(RoadRunnerAPITestsWithLLJit, getIndependentFloatingSpeciesIds){
     getIndependentFloatingSpeciesIds();
 }
diff --git a/test/cxx_api_tests/RoadRunnerAPITestsWithMCJit.cpp b/test/cxx_api_tests/RoadRunnerAPITestsWithMCJit.cpp
index 464fc96e77..546bd09199 100644
--- a/test/cxx_api_tests/RoadRunnerAPITestsWithMCJit.cpp
+++ b/test/cxx_api_tests/RoadRunnerAPITestsWithMCJit.cpp
@@ -194,6 +194,10 @@ TEST_F(RoadRunnerAPITestsWithMCJit, getSelectedValues){
     getSelectedValues();
 }
 
+TEST_F(RoadRunnerAPITestsWithMCJit, getSelectedReactionRates){
+    getSelectedReactionRates();
+}
+
 TEST_F(RoadRunnerAPITestsWithMCJit, getIndependentFloatingSpeciesIds){
     getIndependentFloatingSpeciesIds();
 }