diff --git a/experiment/symbolic_utils.py b/experiment/symbolic_utils.py
index c46c49ed..3d53670b 100644
--- a/experiment/symbolic_utils.py
+++ b/experiment/symbolic_utils.py
@@ -172,17 +172,49 @@ def clean_pred_model(model_str, dataset, est_name):
 
     local_dict = {k: Symbol(k) for k in features}
     new_model_str = model_str
-    # rename features
-    for i, f in reversed(list(enumerate(features))):
-        print("replacing feature", i, "with", f)
-        if any([n in est_name.lower() for n in ["mrgp", "operon", "dsr"]]):
-            i = i + 1
-        new_model_str = new_model_str.replace("x" + str(i), f)
-        new_model_str = new_model_str.replace("x_" + str(i), f)
-        new_model_str = new_model_str.replace("X_" + str(i), f)
-        new_model_str = new_model_str.replace("X" + str(i), f)
-        new_model_str = new_model_str.replace("x[:,{}]".format(i), f)
-        new_model_str = new_model_str.replace("x[{}]".format(i), f)
+    # Rename generic feature placeholders. Model strings are expected to use
+    # dataset feature names, and some datasets genuinely have names such as
+    # ``x1`` and ``x2``. In that case, treating the names as placeholders would
+    # corrupt an already-correct model.
+    placeholder_pattern = r"\b(?:x_|X_|x|X)(\d+)\b"
+    placeholder_tokens = {
+        match.group(0) for match in re.finditer(placeholder_pattern, new_model_str)
+    }
+    feature_names = {str(feature) for feature in features}
+
+    index_offset = int(  # Some methods start counting features at 1
+        any(name in est_name.lower() for name in ["mrgp", "operon", "dsr"])
+    )
+
+    def feature_name(index):
+        feature_index = index - index_offset
+        if 0 <= feature_index < len(features):
+            return str(features[feature_index])
+        return None
+
+    # If every placeholder-looking identifier is an actual feature name, the
+    # string is already in the documented dataset-name form. This intentionally
+    # favors preserving a compliant model in the otherwise ambiguous x1/x2 case.
+    if not placeholder_tokens.issubset(feature_names):
+
+        def replace_placeholder(match):
+            replacement = feature_name(int(match.group(1)))
+            return replacement if replacement is not None else match.group(0)
+
+        new_model_str = re.sub(placeholder_pattern, replace_placeholder, new_model_str)
+
+    # Array-style placeholders cannot be dataset feature names, so they can
+    # always be normalized. Keep this separate from identifier replacement.
+    def replace_array_placeholder(match):
+        replacement = feature_name(int(match.group(1)))
+        return replacement if replacement is not None else match.group(0)
+
+    new_model_str = re.sub(
+        r"x\[:,\s*(\d+)\s*\]", replace_array_placeholder, new_model_str
+    )
+    new_model_str = re.sub(
+        r"x\[\s*(\d+)\s*\]", replace_array_placeholder, new_model_str
+    )
     # operators
     new_model_str = new_model_str.replace("^", "**")
     # GP-GOMEA
diff --git a/experiment/test_algorithm.py b/experiment/test_algorithm.py
index ebae2bf3..ac7c5c74 100644
--- a/experiment/test_algorithm.py
+++ b/experiment/test_algorithm.py
@@ -6,18 +6,20 @@
 from os.path import abspath
 from os.path import dirname as d
 
+import numpy as np
+import pandas as pd
 import pytest
 
 sys.path.append(d(abspath(__file__)))
 
-
+import symbolic_utils
 from evaluate_model import evaluate_model
 from read_file import read_file
 from symbolic_utils import (
     complexity,
     round_floats,
 )
-from sympy import Symbol
+from sympy import Symbol, simplify
 
 # symbolic model stuff
 from sympy.parsing.sympy_parser import parse_expr
@@ -36,6 +38,51 @@
 random_state = 42
 
 
+def _use_toy_dataframe(monkeypatch, feature_names):
+    """Make clean_pred_model see a small dataset without a fixture on disk."""
+    dataframe = pd.DataFrame(np.zeros((2, len(feature_names))), columns=feature_names)
+    dataframe["target"] = 0.0
+
+    def read_toy_file(_dataset):
+        return (
+            dataframe.drop(columns="target"),
+            dataframe["target"].to_numpy(),
+            dataframe.drop(columns="target").columns.to_numpy(),
+        )
+
+    monkeypatch.setattr(symbolic_utils, "read_file", read_toy_file)
+
+
+def test_clean_pred_model_preserves_x_named_features(monkeypatch, ml):
+    feature_names = ["x1", "x2", "x3", "y1", "y2", "y3"]
+    model = "x1*y1 + x2*y2 + x3*y3"
+    _use_toy_dataframe(monkeypatch, feature_names)
+
+    cleaned = symbolic_utils.clean_pred_model(model, "toy.tsv", ml)
+
+    assert simplify(cleaned - parse_expr(model)) == 0
+
+
+def test_clean_pred_model_preserves_mixed_x_named_features(monkeypatch, ml):
+    feature_names = ["m1", "m2", "G", "x1", "x2", "y1", "y2", "z1", "z2"]
+    model = "x1**2 - 2*x1*x2 + x2**2"
+    _use_toy_dataframe(monkeypatch, feature_names)
+
+    cleaned = symbolic_utils.clean_pred_model(model, "toy.tsv", ml)
+
+    assert simplify(cleaned - parse_expr(model)) == 0
+
+
+def test_clean_pred_model_still_normalizes_generic_placeholders(monkeypatch, ml):
+    _use_toy_dataframe(monkeypatch, ["mass", "velocity"])
+
+    is_one_based = any(name in ml.lower() for name in ["mrgp", "operon", "dsr"])
+    model = "x1 + x2" if is_one_based else "x0 + x1"
+    cleaned = symbolic_utils.clean_pred_model(model, "toy.tsv", ml)
+
+    assert simplify(cleaned - parse_expr("mass + velocity")) == 0
+
+
 def get_algorithm(ml):
 
     algorithm = importlib.__import__(