diff --git a/compare.jl b/compare.jl
new file mode 100644
index 00000000..8aa9dd1b
--- /dev/null
+++ b/compare.jl
@@ -0,0 +1,45 @@
+using Revise, StructuralIdentifiability, ParamPunPam, RationalFunctionFields
+import StructuralIdentifiability: parent_ring_change, lie_derivative
+import RationalFunctionFields: fields_equal
+
+cmp = (f, g) -> begin
+    f = parent_ring_change(f, parent(ode))
+    g = parent_ring_change(g, parent(ode))
+    df = lie_derivative(f, ode)
+    dg = lie_derivative(g, ode)
+    default_cmp = RationalFunctionFields.rational_function_cmp
+    if iszero(df) && iszero(dg)
+        default_cmp(f, g)
+    else
+        default_cmp(df, dg)
+    end    
+end
+
+ode = @ODEmodel(
+            x1'(t) =
+                (p1 * x4(t)) - (p3 * x1(t)) -
+                p4 * ((x1(t)^2 / (p5 + x1(t))) * (1 + (p6 * u1(t) / (p7 + u1(t))))),
+            x2'(t) =
+                p8 - (p9 * x2(t)) -
+                p10 * (
+                (x1(t) * x2(t) / (p11 + x2(t))) * (1 + (p12 * u1(t) / (p13 + u1(t))))
+            ),
+            x3'(t) =
+                p14 - (p15 * x3(t)) -
+                p16 * x1(t) * x3(t) * (1 - p18 * u1(t)) / (p17 + x3(t)),
+            x4'(t) =
+                p20 - p21 * (1 - p24) * (1 - p25) / ((p22^4) + 1) - (p20 * x4(t)) +
+                (p21 * (x3(t)^4)) +
+                (1 + p23 * u1(t)) * (1 - p24 * x1(t)) * (1 - p25 * x2(t)) /
+                (p22^4 + x3(t)^4),
+            y1(t) = x1(t),
+            y2(t) = x2(t),
+            y3(t) = x3(t),
+            y4(t) = x4(t)
+        )
+
+id1 = find_identifiable_functions(ode, with_states=true, cmp = cmp)
+
+id2 = find_identifiable_functions(ode, with_states=true, cmp = RationalFunctionFields.rational_function_cmp)
+
+@assert fields_equal(RationalFunctionField(id1), RationalFunctionField(id2), 0.99)
diff --git a/src/global_identifiability.jl b/src/global_identifiability.jl
index 4f97226f..8d0a174e 100644
--- a/src/global_identifiability.jl
+++ b/src/global_identifiability.jl
@@ -91,6 +91,7 @@ The function returns a tuple containing the following:
         with_states::Bool = false,
         var_change_policy = :default,
         rational_interpolator = :VanDerHoevenLecerf,
+        cmp = RationalFunctionFields.rational_function_cmp
     ) where {T}
     @info "Computing IO-equations"
     ioeq_time = @elapsed io_equations =
@@ -148,6 +149,7 @@ The function returns a tuple containing the following:
             seed = 42,
             simplify = :standard,
             rational_interpolator = rational_interpolator,
+            cmp = cmp
         )
         dennums_simplified = fractions_to_dennums(no_states_simplified)
         # switch back the ring
@@ -202,7 +204,7 @@ Output: a list L of booleans with L[i] being the identifiability status of the i
         known = known,
         prob_threshold = half_p,
         var_change_policy = var_change_policy,
-        with_states = states_needed,
+        with_states = states_needed
     )
 
     funcs_to_check = Vector{Generic.FracFieldElem{P}}(
diff --git a/src/identifiable_functions.jl b/src/identifiable_functions.jl
index 02079564..1447123c 100644
--- a/src/identifiable_functions.jl
+++ b/src/identifiable_functions.jl
@@ -54,6 +54,7 @@ function find_identifiable_functions(
         simplify = :standard,
         rational_interpolator = :VanDerHoevenLecerf,
         loglevel = Logging.Info,
+        cmp = RationalFunctionFields.rational_function_cmp,
     ) where {T <: MPolyRingElem{QQFieldElem}}
     restart_logging(loglevel = loglevel)
     reset_timings()
@@ -66,6 +67,7 @@ function find_identifiable_functions(
                 with_states = with_states,
                 simplify = simplify,
                 rational_interpolator = rational_interpolator,
+                cmp = cmp
             )
         else
             id_funcs = _find_identifiable_functions_kic(
@@ -75,6 +77,7 @@ function find_identifiable_functions(
                 seed = seed,
                 simplify = simplify,
                 rational_interpolator = rational_interpolator,
+                cmp = cmp
             )
             # renaming variables from `x(t)` to `x(0)`
             return replace_with_ic(ode, id_funcs)
@@ -89,6 +92,7 @@ function _find_identifiable_functions(
         with_states = false,
         simplify = :standard,
         rational_interpolator = :VanDerHoevenLecerf,
+        cmp = RationalFunctionFields.rational_function_cmp
     ) where {T <: MPolyRingElem{QQFieldElem}}
     Random.seed!(seed)
     @assert simplify in (:standard, :weak, :strong, :absent)
@@ -109,6 +113,7 @@ function _find_identifiable_functions(
         prob_threshold = half_p,
         with_states = with_states,
         rational_interpolator = rational_interpolator,
+        cmp = cmp
     )
     # If simplification is needed
     if simplify !== :absent
@@ -122,7 +127,8 @@ function _find_identifiable_functions(
             seed = seed,
             simplify = simplify,
             rational_interpolator = rational_interpolator,
-            priority_variables = [parent_ring_change(p, bring) for p in ode.parameters],
+            # priority_variables = [parent_ring_change(p, bring) for p in ode.parameters],
+            cmp = cmp
         )
     else
         id_funcs_fracs = RationalFunctionFields.dennums_to_fractions(id_funcs)
diff --git a/src/known_ic.jl b/src/known_ic.jl
index 83f359bc..bd0cccf0 100644
--- a/src/known_ic.jl
+++ b/src/known_ic.jl
@@ -33,6 +33,7 @@ function _find_identifiable_functions_kic(
         seed = 42,
         simplify = :standard,
         rational_interpolator = :VanDerHoevenLecerf,
+        cmp = RationalFunctionFields.rational_function_cmp
     ) where {T <: MPolyRingElem{Nemo.QQFieldElem}}
     Random.seed!(seed)
     @assert simplify in (:standard, :weak, :strong, :absent)
@@ -45,6 +46,7 @@ function _find_identifiable_functions_kic(
         simplify = :absent,
         rational_interpolator = rational_interpolator,
         seed = seed,
+        cmp = cmp
     )
 
     id_funcs = vcat(id_funcs_general, [f // one(parent(ode)) for f in known_ic])
@@ -56,6 +58,7 @@ function _find_identifiable_functions_kic(
             seed = seed,
             simplify = simplify,
             rational_interpolator = rational_interpolator,
+            cmp = cmp
         )
     end