diff --git a/share/lib/python/neuron/rxd/generalizedReaction.py b/share/lib/python/neuron/rxd/generalizedReaction.py
index e20a6c9274..a7914c95ad 100644
--- a/share/lib/python/neuron/rxd/generalizedReaction.py
+++ b/share/lib/python/neuron/rxd/generalizedReaction.py
@@ -44,7 +44,8 @@ def get_scheme_rate1_rate2_regions_custom_dynamics_mass_action(args, kwargs):
             rate1 = args[1]
             rate2 = args[2]
         else:
-            scheme = args[0] > args[1]
+            scheme = args[0].__gt__(args[1])
+            # don't use '>' to avoid Python prioritizing subclass methods
             rate1 = args[2]
             rate2 = None
     elif len(args) == 2:
diff --git a/share/lib/python/neuron/rxd/reaction.py b/share/lib/python/neuron/rxd/reaction.py
index d98206c9ce..f5c9d9aff3 100644
--- a/share/lib/python/neuron/rxd/reaction.py
+++ b/share/lib/python/neuron/rxd/reaction.py
@@ -102,27 +102,18 @@ def _do_init(self):
     def _update_rates(self):
         lhs = self._scheme._lhs._items
         rhs = self._scheme._rhs._items
-        if self._dir == "<":
-            # TODO: remove this limitation (probably means doing with rate_b what done with rate_f and making sure _sources and _dests are correct
-            raise RxDException(
-                "pure reverse reaction currently not supported; reformulate as a forward reaction"
-            )
 
         rate_f = rxdmath._ensure_arithmeticed(self._original_rate_f)
         rate_b = rxdmath._ensure_arithmeticed(self._original_rate_b)
 
         if not self._custom_dynamics:
-            for k, v in zip(list(lhs.keys()), list(lhs.values())):
-                if v == 1:
-                    rate_f *= k
-                else:
-                    rate_f *= k**v
-            if self._dir == "<>":
+            if ">" in self._dir:
+                for k, v in lhs.items():
+                    rate_f *= k if v == 1 else k**v
+            if "<" in self._dir:
                 for k, v in rhs.items():
-                    if v == 1:
-                        rate_b *= k
-                    else:
-                        rate_b *= k**v
+                    rate_b *= k if v == 1 else k**v
+
         rate = rate_f - rate_b
         self._rate_arithmeticed = rate
 
diff --git a/test/rxd/test_state_reaction.py b/test/rxd/test_state_reaction.py
new file mode 100644
index 0000000000..c3b09dd5fb
--- /dev/null
+++ b/test/rxd/test_state_reaction.py
@@ -0,0 +1,19 @@
+import pytest
+from .testutils import compare_data, tol
+
+
+def test_state_reaction(neuron_instance):
+    """Test a simple rxd.Reaction involving a rxd.State"""
+    h, rxd, data, save_path = neuron_instance
+    sec = h.Section(name="sec")
+    cyt = rxd.Region([sec], name="cyt", nrn_region="i")
+    ca = rxd.Species(cyt, name="ca", charge=2, initial=0.001)
+    cabuf = rxd.State(cyt, name="cabuf", initial=0)
+    # buffering = rxd.Reaction(ca, cabuf, 0.02)
+    # generated the test data with a forward reaction
+    buffering = rxd.Reaction(cabuf < ca, 0.02)
+    h.finitialize(-65)
+    h.continuerun(10)
+    if not save_path:
+        max_err = compare_data(data)
+        assert max_err < tol
diff --git a/test/rxd/testdata b/test/rxd/testdata
index 1512df7c1b..002556b624 160000
--- a/test/rxd/testdata
+++ b/test/rxd/testdata
@@ -1 +1 @@
-Subproject commit 1512df7c1bab225f9abc2f30cfd19c084d814d1e
+Subproject commit 002556b62494846fa3788e6bb87a8f519c60859e