diff --git a/devito/data/utils.py b/devito/data/utils.py
index d221978eaa..45f6c376bc 100644
--- a/devito/data/utils.py
+++ b/devito/data/utils.py
@@ -47,18 +47,20 @@ def index_apply_modulo(idx, modulo):
     if is_integer(idx):
         return idx % modulo
     elif isinstance(idx, slice):
+        # A negative bound counts from the end, so it is its magnitude that
+        # wraps around; e.g., `-1:` must select the last entry whatever `modulo`
         if idx.start is None:
             start = idx.start
         elif idx.start >= 0:
             start = idx.start % modulo
         else:
-            start = -(idx.start % modulo)
+            start = -(-idx.start % modulo)
         if idx.stop is None:
             stop = idx.stop
         elif idx.stop >= 0:
             stop = idx.stop % (modulo + 1)
         else:
-            stop = -(idx.stop % (modulo + 1))
+            stop = -(-idx.stop % (modulo + 1))
         return slice(start, stop, idx.step)
     elif isinstance(idx, (tuple, list)):
         return [i % modulo for i in idx]
diff --git a/tests/test_data.py b/tests/test_data.py
index 047af71c0c..28698b9344 100644
--- a/tests/test_data.py
+++ b/tests/test_data.py
@@ -93,6 +93,28 @@ def test_negative_start(self):
         f.data[idx] = dat
         assert np.all(np.array(f.data[9:]) == dat)
 
+    @pytest.mark.parametrize('idx', [
+        slice(-1, None), slice(-2, None), slice(-3, None), slice(None, -1),
+        slice(None, -2), slice(1, -1), slice(-2, -1), slice(-1, None, -1)
+    ])
+    def test_negative_slicing_modulo(self, idx):
+        """
+        Test slicing with negative bounds along a modulo-indexed Dimension, here
+        the time Dimension of a buffered TimeFunction. Both reads and writes must
+        select the same entries as NumPy does.
+        """
+        grid = Grid(shape=(4,))
+        u = TimeFunction(name='u', grid=grid, time_order=2)
+        assert u.data._modulo == (True, False)
+
+        expected = np.arange(12, dtype=np.float32).reshape(3, 4)
+        u.data[:] = expected
+        assert np.array_equal(u.data[idx], expected[idx])
+
+        expected[idx] = -1.
+        u.data[idx] = -1.
+        assert np.array_equal(u.data, expected)
+
     def test_halo_indexing(self):
         """Test data packing/unpacking in presence of a halo region."""
         domain_shape = (16, 16, 16)