diff --git a/test/unit/test_memory.py b/test/unit/test_memory.py
index 1dade5bf3..5f54320a8 100644
--- a/test/unit/test_memory.py
+++ b/test/unit/test_memory.py
@@ -6,6 +6,16 @@
     pauli_term_to_measurement_memory_map,
 )
 from pyquil.paulis import sX, sY
+from pyquil.quilatom import (
+    MemoryReference,
+    quil_cis,
+    quil_cos,
+    quil_exp,
+    quil_sin,
+    quil_sqrt,
+    substitute,
+    substitute_array,
+)
 
 
 def test_merge_memory_map_lists():
@@ -34,3 +44,42 @@ def test_pauli_term_to_measurement_memory_map():
         "measurement_beta": [0.0, np.pi / 2],
         "measurement_gamma": [0.0, -np.pi / 2],
     }
+
+
+def test_substitute_memory_reference():
+    x = MemoryReference("x", 0, declared_size=2)
+    x_0 = MemoryReference("x", 0, declared_size=2)
+    x_1 = MemoryReference("x", 1, declared_size=2)
+
+    # complete substitutions
+
+    assert substitute(x_0, {x: [+5, -5]}) == +5
+    assert substitute(x_1, {x: [+5, -5]}) == -5
+
+    assert substitute(x_0 + x_1, {x: [+5, -5]}) == 0
+    assert substitute(x_0 - x_1, {x: [+5, -5]}) == 10
+    assert substitute(x_0 * x_1, {x: [+5, -5]}) == -25
+    assert substitute(x_0 / x_1, {x: [+5, -5]}) == -1
+
+    assert substitute(x_0 * x_0 ** 2 / x_1, {x: [5, 10]}) == 12.5
+
+    assert np.isclose(substitute(quil_exp(x_0), {x: [5, 10]}), np.exp(5))
+    assert np.isclose(substitute(quil_sin(x_0), {x: [5, 10]}), np.sin(5))
+    assert np.isclose(substitute(quil_cos(x_0), {x: [5, 10]}), np.cos(5))
+    assert np.isclose(substitute(quil_sqrt(x), {x: [5, 10]}), np.sqrt(5))
+    assert np.isclose(substitute(quil_cis(x), {x: [5, 10]}), np.exp(1j * 5.0))
+
+    # incomplete substitutions
+
+    y = MemoryReference("y", 0, declared_size=1)
+    z = MemoryReference("z", 0, declared_size=1)
+
+    assert substitute(y + z, {y: [5]}) == 5 + z
+
+    assert substitute(quil_cis(z), {y: [5]}) == quil_cis(z)
+
+    # array substitution pass-through
+
+    a = MemoryReference("a", 0, declared_size=1)
+
+    assert np.allclose(substitute_array([quil_sin(a), quil_cos(a)], {a: [5]}), [np.sin(5), np.cos(5)])