diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index dd2e34c..aef4240 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -25,7 +25,7 @@ jobs:
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: julia-actions/cache@v1
+      - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@latest
       - run: |
           git config --global user.name Tester
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index 1b195d8..7c06bf8 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -19,7 +19,7 @@ jobs:
       - uses: julia-actions/setup-julia@v2
         with:
           version: '1.10'
-      - uses: julia-actions/cache@v1
+      - uses: julia-actions/cache@v2
       - name: Install dependencies
         run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
       - name: Build and deploy
diff --git a/.gitignore b/.gitignore
index 7b355c3..dfc63f2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,3 @@
-/docs/build/
 
 # ignore vscode
 ./vscode
@@ -12,9 +11,7 @@
 /docs/build
 Manifest.toml
 
-dev/Alloascoidea_hylecoeti.json
-dev/Alloascoidea_hylecoeti_isozymes.json
-dev/Alloascoidea_hylecoeti_molar_mass.json
-
-
 *.json
+
+# vim
+.*.swp
diff --git a/Project.toml b/Project.toml
index 674c611..701cf8d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "DifferentiableMetabolism"
 uuid = "9677e053-0fc1-49e6-98a8-2ad3a1525e54"
 authors = ["St. Elmo Wilken <stelmozors@gmail.com> and contributors"]
-version = "0.2.0"
+version = "0.3.0"
 
 [deps]
 AbstractFBCModels = "5a4f3dfa-1789-40f8-8221-69268c29937c"
diff --git a/dev/dev.jl b/dev/dev.jl
deleted file mode 100644
index 91ceb6e..0000000
--- a/dev/dev.jl
+++ /dev/null
@@ -1,105 +0,0 @@
-using COBREXA, JSON, HiGHS, FastDifferentiation, DifferentiableMetabolism
-import AbstractFBCModels.CanonicalModel as CM
-import AbstractFBCModels as A
-import JSONFBCModels
-
-
-mass_bound = 350.0
-optimizer = HiGHS.Optimizer
-gene_zero_tol = 1e-8
-flux_zero_tol = 1e-8
-
-
-model = convert(CM.Model, load_model("dev/Alloascoidea_hylecoeti.json"));
-reaction_isozymes = Dict(
-    k => Dict(
-        "$(k)_$i" => Isozyme(
-            d["gene_product_stoichiometry"],
-            d["kcat_forward"],
-            d["kcat_reverse"],
-        ) for (i, d) in enumerate(v)
-    ) for (k, v) in JSON.parsefile("dev/Alloascoidea_hylecoeti_isozymes.json")
-);
-gene_product_molar_masses =
-    Dict(x => y for (x, y) in JSON.parsefile("dev/Alloascoidea_hylecoeti_molar_mass.json"))
-capacity = [("uncategorized", A.genes(model), mass_bound)]
-# open exchange reactions
-for (r, rxn) in model.reactions
-    if contains(rxn.name, "exchange")
-        model.reactions[r].lower_bound = rxn.lower_bound * 1000.0
-        model.reactions[r].upper_bound = rxn.upper_bound * 1000.0
-    end
-end
-
-ec_solution = enzyme_constrained_flux_balance_analysis(
-    model;
-    reaction_isozymes,
-    gene_product_molar_masses,
-    capacity,
-    optimizer,
-)
-
-reaction_parameter_isozymes = Dict(
-    r => Dict(
-        k => ParameterIsozyme(
-            v.gene_product_stoichiometry,
-            FastDifferentiation.Node(Symbol("$(k)_f")),
-            FastDifferentiation.Node(Symbol("$(k)_r")),
-        ) for (k, v) in isos
-    ) for (r, isos) in reaction_isozymes
-);
-
-pruned_reaction_isozymes =
-    prune_reaction_isozymes(reaction_parameter_isozymes, ec_solution, gene_zero_tol);
-
-pruned_gene_product_molar_masses =
-    prune_gene_product_molar_masses(gene_product_molar_masses, ec_solution, gene_zero_tol);
-
-pruned_model = prune_model(model, ec_solution, flux_zero_tol, gene_zero_tol)
-#make reactions uni-directional
-for (r,rxn) in pruned_model.reactions
-    if ec_solution.fluxes[r] < 0 
-        pruned_model.reactions[r].upper_bound = 0.0
-    else
-        pruned_model.reactions[r].lower_bound = 0.0
-    end
-end
-
-
-pkm = build_kinetic_model(
-    pruned_model;
-    reaction_isozymes=pruned_reaction_isozymes,
-    gene_product_molar_masses=pruned_gene_product_molar_masses,
-    capacity,
-);
-
-parameter_values = Dict{Symbol,Float64}();
-for (r, iso) in pruned_reaction_isozymes
-    for (k, v) in iso
-        parameter_values[Symbol(k, :_f)] = reaction_isozymes[r][k].kcat_forward
-        parameter_values[Symbol(k, :_r)] = reaction_isozymes[r][k].kcat_reverse
-    end
-end
-parameters = collect(keys(parameter_values));
-
-# costly preparation of KKTs
-pkmKKT, vids = differentiate_prepare_kkt(pkm, pkm.objective.value, parameters)
-
-ec_solution2, x_vals, eq_dual_vals, ineq_dual_vals = optimized_constraints_with_parameters(
-    pkm,
-    parameter_values;
-    objective=pkm.objective.value,
-    optimizer=HiGHS.Optimizer,
-    modifications=[silence],
-);
-
-# cheap symbolic differentiation
-sens = differentiate_solution(
-    pkmKKT,
-    x_vals,
-    eq_dual_vals,
-    ineq_dual_vals,
-    parameter_values,
-    parameters;
-    scale=true, # unitless sensitivities
-);
diff --git a/docs/Project.toml b/docs/Project.toml
index 19ec310..134fdd4 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -7,10 +7,10 @@ ConstraintTrees = "5515826b-29c3-47a5-8849-8513ac836620"
 DifferentiableMetabolism = "9677e053-0fc1-49e6-98a8-2ad3a1525e54"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+FastDifferentiation = "eb9bf01b-bf85-4b60-bf87-ee5de06c00be"
 JSONFBCModels = "475c1105-d6ed-49c1-9b32-c11adca6d3e8"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
-Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Tulip = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa"
 
 [compat]
diff --git a/docs/src/1-parametric-models.jl b/docs/src/1-parametric-models.jl
index bb7160a..c8007ce 100644
--- a/docs/src/1-parametric-models.jl
+++ b/docs/src/1-parametric-models.jl
@@ -16,7 +16,7 @@
 
 using DifferentiableMetabolism
 
-using Symbolics
+using FastDifferentiation
 using ConstraintTrees
 using COBREXA
 using Tulip
@@ -46,7 +46,7 @@ base_model.fluxes
 # ## Add parameters to the model
 
 # Make bound of r2 and mass balance of m3 parameters
-Symbolics.@variables r2bound m3bound
+@variables r2bound m3bound
 
 m.fluxes.r2 =
     ConstraintTrees.Constraint(m.fluxes.r2.value, -2 * ParameterBetween(r2bound, 0))
@@ -55,7 +55,7 @@ m.flux_stoichiometry.m3 =
     ConstraintTrees.Constraint(m.flux_stoichiometry.m3.value, ParameterEqualTo(m3bound) / 2)
 
 # # add parametric constraints
-Symbolics.@variables p[1:4]
+p = make_variables(:p, 4)
 
 m *=
     :linparam^ConstraintTrees.Constraint(
@@ -65,11 +65,11 @@ m *=
 
 # substitute params into model
 parameter_substitutions = Dict(
-    r2bound => 4.0,
-    m3bound => 0.1, # lose some mass here
-    p[1] => 1.0,
-    p[2] => 1.0,
-    p[3] => 4.0,
+    :r2bound => 4.0,
+    :m3bound => 0.1, # lose some mass here
+    :p1 => 1.0,
+    :p2 => 1.0,
+    :p3 => 4.0,
 )
 
 m_noparams, _, _, _ = optimized_constraints_with_parameters(
@@ -85,7 +85,7 @@ m_noparams.fluxes
 # ## Change the parameters and re-solve
 
 # substitute params into model
-parameter_substitutions[m3bound] = 0.0
+parameter_substitutions[:m3bound] = 0.0
 
 m_noparams, _, _, _ = optimized_constraints_with_parameters(
     m,
@@ -99,7 +99,7 @@ m_noparams.fluxes
 
 # ## Quadratic parameters also work
 
-Symbolics.@variables q[1:6]
+q = make_variables(:q, 6)
 
 m.objective = ConstraintTrees.Constraint(
     value = sum(
@@ -110,7 +110,8 @@ m.objective = ConstraintTrees.Constraint(
 
 m *= :objective_bound^ConstraintTrees.Constraint(value = m.fluxes.r6.value, bound = 2.0)
 
-parameter_substitutions = merge(parameter_substitutions, Dict(zip(q, fill(1.0, 6))))
+parameter_substitutions =
+    merge(parameter_substitutions, Dict(v.node_value => 1.0 for v in q))
 
 m_noparams, _, _, _ = optimized_constraints_with_parameters(
     m,
@@ -121,138 +122,136 @@ m_noparams, _, _, _ = optimized_constraints_with_parameters(
 )
 m_noparams.fluxes
 
+v(x) = DifferentiableMetabolism.substitute(x, _ -> throw("no subst!")) #src
 @test isapprox(m_noparams.objective, 11.0; atol = TEST_TOLERANCE) #src
 pb = ParameterBetween(1, 2.0) #src
-@test (-pb).lower == -1 && (-pb).upper == -2  #src
-@test (pb * 2).lower == 2 && (pb * 2).upper == 4  #src
-@test (pb / 2).lower == 0.5 && (pb / 2).upper == 1  #src
+@test v((-pb).lower) == -1 && v((-pb).upper) == -2 #src
+@test v((pb * 2).lower) == 2 && v((pb * 2).upper) == 4 #src
+@test v((pb / 2).lower) == 0.5 && v((pb / 2).upper) == 1 #src
 pe = ParameterEqualTo(1) #src
-@test (-pe).equal_to == -1 #src
-@test (pe * 2).equal_to == 2 #src
-@test (pe / 2).equal_to == 0.5 #src
+@test v((-pe).equal_to) == -1 #src
+@test v((pe * 2).equal_to) == 2 #src
+@test v((pe / 2).equal_to) == 0.5 #src
 plv1 = ParameterLinearValue(ConstraintTrees.LinearValue(1)) #src
-@test all(plv1.idxs .== [0]) && all(plv1.weights .== [Num(1.0)]) #src
+@test plv1.idxs == [0] && all(v.(plv1.weights) .== 1.0) #src
 plv2 = ParameterLinearValue([0], [1.0]) #src
-@test all(plv2.idxs .== [0]) && all(plv2.weights .== [Num(1.0)]) #src
+@test plv2.idxs == [0] && all(v.(plv2.weights) .== 1.0) #src
 plv3 = ParameterLinearValue(0) #src
 @test isempty(plv3.idxs) && isempty(plv3.weights) #src
 plv4 = ParameterLinearValue(1) #src
-@test all(plv4.idxs .== [0]) && all(plv4.weights .== [Num(1.0)]) #src
+@test plv4.idxs == [0] && all(v.(plv4.weights) .== 1.0) #src
 plv5 = convert(ParameterLinearValue, 1) #src
-@test all(plv5.idxs .== [0]) && all(plv5.weights .== [Num(1.0)]) #src
+@test plv5.idxs == [0] && all(v.(plv5.weights) .== 1.0) #src
 plv6 = zero(ParameterLinearValue) #src
 @test isempty(plv6.idxs) && isempty(plv6.weights) #src
 plv7 = plv5 + 1 #src
-@test all(plv7.idxs .== [0]) && all(plv7.weights .== [Num(2.0)]) #src
+@test all(v.(plv7.idxs) .== 0) && all(v.(plv7.weights) .== 2.0) #src
 plv8 = ParameterLinearValue([1, 2, 3], [1, 2, 3]) #src
 plv9 = ParameterLinearValue([1, 3], [4.0, 5.0]) #src
 plv10 = plv8 + plv9 #src
-@test all(plv10.idxs .== [1, 2, 3]) && all(plv10.weights .== [Num(5), Num(2), Num(8)]) #src
+@test plv10.idxs == [1, 2, 3] && v.(plv10.weights) == [5, 2, 8] #src
 plv10 = plv9 + plv8 #src
-@test all(plv10.idxs .== [1, 2, 3]) && all(plv10.weights .== [Num(5), Num(2), Num(8)]) #src
+@test plv10.idxs == [1, 2, 3] && v.(plv10.weights) == [5, 2, 8] #src
 plv11 = plv5 - 2 #src
-@test all(plv11.idxs .== [0]) && all(plv11.weights .== [Num(-1)]) #src
+@test plv11.idxs == [0] && v.(plv11.weights) == [-1] #src
 plv12 = 2 - plv5 #src
-@test all(plv12.idxs .== [0]) && all(plv12.weights .== [Num(1)]) #src
+@test plv12.idxs == [0] && v.(plv12.weights) == [1] #src
 plv13 = plv5 / -2 #src
-@test all(plv13.idxs .== [0]) && all(plv13.weights .== [Num(-0.5)]) #src
+@test plv13.idxs == [0] && v.(plv13.weights) == [-0.5] #src
 plv14 = plv8 - plv9 #src
-@test all(plv14.idxs .== [1, 2, 3]) && all(plv14.weights .== [Num(-3), Num(2), Num(-2)]) #src
+@test plv14.idxs == [1, 2, 3] && v.(plv14.weights) == [-3, 2, -2] #src
 plv15 = 1.0 + ParameterLinearValue(0) #src
-@test all(plv15.idxs .== [0]) && all(plv15.weights .== [Num(1)]) #src
+@test plv15.idxs == [0] && v.(plv15.weights) == [1] #src
 pqv1 = ParameterQuadraticValue([(1, 1)], [1]) #src
 pqv2 = ParameterQuadraticValue(ConstraintTrees.QuadraticValue([(1, 1)], [1])) #src
-@test all(pqv2.idxs[1] .== (1, 1)) && all(pqv2.weights .== Num(1)) #src
+@test pqv2.idxs == [(1, 1)] && v.(pqv2.weights) == [1] #src
 pqv3 = ParameterQuadraticValue(0) #src
 @test isempty(pqv3.idxs) && isempty(pqv3.weights) #src
 pqv4 = ParameterQuadraticValue(ParameterLinearValue([1], [1])) #src
-@test all(pqv4.idxs .== [(0, 1)]) && all(pqv4.weights .== [Num(1)]) #src
+@test pqv4.idxs == [(0, 1)] && v.(pqv4.weights) == [1] #src
 pqv5 = ParameterQuadraticValue(1) #src
-@test all(pqv5.idxs .== [(0, 0)]) && all(pqv5.weights .== [Num(1)]) #src
+@test pqv5.idxs == [(0, 0)] && v.(pqv5.weights) == [1] #src
 pqv6 = convert(ParameterQuadraticValue, 1) #src
-@test all(pqv6.idxs .== [(0, 0)]) && all(pqv6.weights .== [Num(1)]) #src
+@test pqv6.idxs == [(0, 0)] && v.(pqv6.weights) == [1] #src
 pqv7 = convert(ParameterQuadraticValue, ParameterLinearValue([1], [1])) #src
-@test all(pqv7.idxs .== [(0, 1)]) && all(pqv7.weights .== [Num(1)]) #src
+@test pqv7.idxs == [(0, 1)] && v.(pqv7.weights) == [1] #src
 pqv8 = zero(ParameterQuadraticValue) #src
 @test isempty(pqv8.idxs) && isempty(pqv8.weights) #src
 pqv9 = 1 + pqv1 #src
-@test all(pqv9.idxs .== [(0, 0), (1, 1)]) && all(pqv9.weights .== [Num(1), Num(1)]) #src
+@test pqv9.idxs == [(0, 0), (1, 1)] && v.(pqv9.weights) == [1, 1] #src
 pqv10 = ParameterLinearValue([1], [1]) + pqv1 #src
-@test all(pqv10.idxs .== [(0, 1), (1, 1)]) && all(pqv10.weights .== [Num(1), Num(1)]) #src
+@test pqv10.idxs == [(0, 1), (1, 1)] && v.(pqv10.weights) == [1, 1] #src
 pqv11 = 3 - pqv1 #src
-@test all(pqv11.idxs .== [(0, 0), (1, 1)]) && all(pqv11.weights .== [Num(3), Num(-1)]) #src
+@test pqv11.idxs == [(0, 0), (1, 1)] && v.(pqv11.weights) == [3, -1] #src
 pqv12 = ParameterLinearValue([1], [1]) - pqv1 #src
-@test all(pqv12.idxs .== [(0, 1), (1, 1)]) && all(pqv12.weights .== [Num(1), Num(-1)]) #src
+@test pqv12.idxs == [(0, 1), (1, 1)] && v.(pqv12.weights) == [1, -1] #src
 pqv13 = pqv1 - ParameterLinearValue([1], [1]) #src
-@test all(pqv13.idxs .== [(0, 1), (1, 1)]) && all(pqv13.weights .== [Num(-1), Num(1)]) #src
+@test pqv13.idxs == [(0, 1), (1, 1)] && v.(pqv13.weights) == [-1, 1] #src
 pqv14 = 42 * pqv1 #src
-@test all(pqv14.idxs .== [(1, 1)]) && all(pqv14.weights .== [Num(42)]) #src
+@test pqv14.idxs == [(1, 1)] && v.(pqv14.weights) == [42] #src
 pqv15 = pqv1 - (3 * pqv1) #src
-@test all(pqv15.idxs .== [(1, 1)]) && all(pqv15.weights .== [Num(-2)]) #src
+@test pqv15.idxs == [(1, 1)] && v.(pqv15.weights) == [-2] #src
 pqv16 = pqv1 / 2 #src
-@test all(pqv16.idxs .== [(1, 1)]) && all(pqv16.weights .== [Num(0.5)]) #src
-pqv17 = ParameterLinearValue([1, 2], [2, 1]) * ParameterLinearValue([2], [1]) #src 
-@test all(pqv17.idxs .== [(1, 2), (2, 2)]) && all(pqv17.weights .== [Num(2), Num(1)]) #src
+@test pqv16.idxs == [(1, 1)] && v.(pqv16.weights) == [0.5] #src
+pqv17 = ParameterLinearValue([1, 2], [2, 1]) * ParameterLinearValue([2], [1]) #src
+@test pqv17.idxs == [(1, 2), (2, 2)] && v.(pqv17.weights) == [2, 1] #src
 pqv18 = ParameterQuadraticValue([(1, 1), (2, 2)], [1.0, 3.0]) #src
 pqv19 = pqv18 + pqv1 #src
-@test all(pqv19.idxs .== [(1, 1), (2, 2)]) && all(pqv19.weights .== [Num(2), Num(3)]) #src
+@test pqv19.idxs == [(1, 1), (2, 2)] && v.(pqv19.weights) == [2, 3] #src
 pqv20 = pqv1 + pqv18 #src
-@test all(pqv20.idxs .== [(1, 1), (2, 2)]) && all(pqv19.weights .== [Num(2), Num(3)]) #src
+@test pqv20.idxs == [(1, 1), (2, 2)] && v.(pqv19.weights) == [2, 3] #src
 pqv21 = ParameterQuadraticValue([(1, 1), (2, 2)], [1, 2]) #src
-@test all(pqv21.idxs .== [(1, 1), (2, 2)]) && all(pqv21.weights .== [Num(1), Num(2)]) #src
+@test pqv21.idxs == [(1, 1), (2, 2)] && v.(pqv21.weights) == [1, 2] #src
 pqv22 = convert(ParameterQuadraticValue, ConstraintTrees.LinearValue([1], [1])) #src
-@test all(pqv22.idxs .== [(0, 1)]) && all(pqv22.weights .== [Num(1)]) #src
+@test pqv22.idxs == [(0, 1)] && v.(pqv22.weights) == [1] #src
 pqv23 = convert( #src
     ParameterQuadraticValue, #src
     ConstraintTrees.QuadraticValue([(1, 1), (2, 2)], [1, 2]), #src
 ) #src
-@test all(pqv23.idxs .== [(1, 1), (2, 2)]) && all(pqv23.weights .== [Num(1), Num(2)]) #src
-pqv24 = Num(5) - ConstraintTrees.QuadraticValue([(1, 1)], [1]) #src
-@test all(pqv24.idxs .== [(0, 0), (1, 1)]) && all(pqv24.weights .== [Num(5), Num(-1)]) #src
+@test pqv23.idxs == [(1, 1), (2, 2)] && v.(pqv23.weights) == [1, 2] #src
+pqv24 = 5 - ConstraintTrees.QuadraticValue([(1, 1)], [1]) #src
+@test pqv24.idxs == [(0, 0), (1, 1)] && v.(pqv24.weights) == [5, -1] #src
 pqv25 = pqv18 - 5.0 #src
-@test all(pqv25.idxs .== [(0, 0), (1, 1), (2, 2)]) && #src
-      all(pqv25.weights .== [Num(-5), Num(1), Num(3)]) #src
-pr1 = Num(1) + ConstraintTrees.LinearValue([1, 2], [1, 2]) #src
-@test all(pr1.idxs .== [0, 1, 2]) && all(pr1.weights .== [Num(1), Num(1), Num(2)]) #src
-pr2 = ConstraintTrees.LinearValue([1, 2], [1, 2]) - Num(3) #src
-@test all(pr2.idxs .== [0, 1, 2]) && all(pr2.weights .== [Num(-3), Num(1), Num(2)]) #src
-pr2b = Num(3) - ConstraintTrees.LinearValue([1, 2], [1, 2]) #src
-@test all(pr2b.idxs .== [0, 1, 2]) && all(pr2b.weights .== [Num(3), Num(-1), Num(-2)]) #src
-pr3 = Num(3) * ConstraintTrees.LinearValue([1, 2], [1, 2]) #src
-@test all(pr3.idxs .== [1, 2]) && all(pr3.weights .== [Num(3), Num(6)]) #src
-pr3 = ConstraintTrees.LinearValue([1, 2], [1, 2]) / Num(2) #src
-@test all(pr3.idxs .== [1, 2]) && all(pr3.weights .== [Num(0.5), Num(1)]) #src
+@test pqv25.idxs == [(0, 0), (1, 1), (2, 2)] && v.(pqv25.weights) == [-5, 1, 3] #src
+pr1 = 1 + ConstraintTrees.LinearValue([1, 2], [1, 2]) #src
+@test pr1.idxs == [0, 1, 2] && v.(pr1.weights) == [1, 1, 2] #src
+pr2 = ConstraintTrees.LinearValue([1, 2], [1, 2]) - 3 #src
+@test pr2.idxs == [0, 1, 2] && v.(pr2.weights) == [-3, 1, 2] #src
+pr2b = 3 - ConstraintTrees.LinearValue([1, 2], [1, 2]) #src
+@test pr2b.idxs == [0, 1, 2] && v.(pr2b.weights) == [3, -1, -2] #src
+pr3 = 3 * ConstraintTrees.LinearValue([1, 2], [1, 2]) #src
+@test pr3.idxs == [1, 2] && v.(pr3.weights) == [3, 6] #src
+pr3 = ConstraintTrees.LinearValue([1, 2], [1, 2]) / 2 #src
+@test pr3.idxs == [1, 2] && v.(pr3.weights) == [0.5, 1] #src
 pr4 = ConstraintTrees.LinearValue([1, 2], [1, 2]) + ParameterLinearValue([1, 2], [1, 2]) #src
-@test all(pr4.idxs .== [1, 2]) && all(pr4.weights .== [Num(2), Num(4)]) #src
+@test pr4.idxs == [1, 2] && v.(pr4.weights) == [2, 4] #src
 pr5 = ParameterLinearValue([1, 2], [2, 3]) - ConstraintTrees.LinearValue([1, 2], [1, 2]) #src
-@test all(pr5.idxs .== [1, 2]) && all(pr5.weights .== [Num(1), Num(1)]) #src
-pr6 = Num(1) + ConstraintTrees.QuadraticValue([(1, 1)], [1]) #src
-@test all(pr6.idxs .== [(0, 0), (1, 1)]) && all(pr6.weights .== [Num(1), Num(1)]) #src
-pr7 = ConstraintTrees.QuadraticValue([(1, 1)], [1]) - Num(1) #src
-@test all(pr7.idxs .== [(0, 0), (1, 1)]) && all(pr7.weights .== [Num(-1), Num(1)]) #src
-pr8 = Num(2) * ConstraintTrees.QuadraticValue([(1, 1)], [1]) #src
-@test all(pr8.idxs .== [(1, 1)]) && all(pr8.weights .== [Num(2)]) #src
-pr9 = ConstraintTrees.QuadraticValue([(1, 1)], [1]) / Num(2)  #src
-@test all(pr9.idxs .== [(1, 1)]) && all(pr9.weights .== [Num(0.5)]) #src
+@test pr5.idxs == [1, 2] && v.(pr5.weights) == [1, 1] #src
+pr6 = 1 + ConstraintTrees.QuadraticValue([(1, 1)], [1]) #src
+@test pr6.idxs == [(0, 0), (1, 1)] && v.(pr6.weights) == [1, 1] #src
+pr7 = ConstraintTrees.QuadraticValue([(1, 1)], [1]) - 1 #src
+@test pr7.idxs == [(0, 0), (1, 1)] && v.(pr7.weights) == [-1, 1] #src
+pr8 = 2 * ConstraintTrees.QuadraticValue([(1, 1)], [1]) #src
+@test pr8.idxs == [(1, 1)] && v.(pr8.weights) == [2] #src
+pr9 = ConstraintTrees.QuadraticValue([(1, 1)], [1]) / 2 #src
+@test pr9.idxs == [(1, 1)] && v.(pr9.weights) == [0.5] #src
 pr10 = #src
     ConstraintTrees.QuadraticValue([(1, 1)], [1]) + ParameterQuadraticValue([(1, 1)], [1]) #src
-@test all(pr10.idxs .== [(1, 1)]) && all(pr10.weights .== [Num(2)]) #src
+@test pr10.idxs == [(1, 1)] && v.(pr10.weights) == [2] #src
 pr11 = #src
     ParameterQuadraticValue([(1, 1)], [2]) - ConstraintTrees.QuadraticValue([(1, 1)], [1]) #src
-@test all(pr11.idxs .== [(1, 1)]) && all(pr11.weights .== [Num(1)]) #src
+@test pr11.idxs == [(1, 1)] && v.(pr11.weights) == [1] #src
 pr12 = ConstraintTrees.LinearValue([1, 2], [1, 2]) + ParameterQuadraticValue([(1, 1)], [2]) #src
-@test all(pr12.idxs .== [(0, 1), (1, 1), (0, 2)]) && #src
-      all(pr12.weights .== [Num(1), Num(2), Num(2)]) #src
+@test pr12.idxs == [(0, 1), (1, 1), (0, 2)] && v.(pr12.weights) == [1, 2, 2] #src
 pr13 = ParameterQuadraticValue([(1, 1)], [2]) - ConstraintTrees.LinearValue([1, 2], [1, 2]) #src
-@test all(pr13.idxs .== [(0, 1), (1, 1), (0, 2)]) && #src
-      all(pr13.weights .== [Num(-1), Num(2), Num(-2)]) #src
+@test pr13.idxs == [(0, 1), (1, 1), (0, 2)] && v.(pr13.weights) == [-1, 2, -2] #src
 pr14 = ParameterLinearValue([1, 2], [2, 3]) * ConstraintTrees.LinearValue([1, 2], [1, 2]) #src
-@test all(pr14.idxs .== [(1, 1), (1, 2), (2, 2)]) && #src
-      all(pr14.weights .== [Num(2), Num(7), Num(6)]) #src
+@test pr14.idxs == [(1, 1), (1, 2), (2, 2)] && v.(pr14.weights) == [2, 7, 6] #src
 pr15 = ConstraintTrees.LinearValue([1, 2], [1, 2]) - ParameterQuadraticValue([(1, 1)], [2]) #src
-@test all(pr15.idxs .== [(0, 1), (1, 1), (0, 2)]) && #src
-      all(pr15.weights .== [Num(1), Num(-2), Num(2)]) #src
-@test ConstraintTrees.substitute( #src
-    ParameterQuadraticValue([(0, 0), (1, 1)], [1, 2]), #src
-    [Num(1), Num(2)], #src
+@test pr15.idxs == [(0, 1), (1, 1), (0, 2)] && v.(pr15.weights) == [1, -2, 2] #src
+@test v( #src
+    ConstraintTrees.substitute( #src
+        ParameterQuadraticValue([(0, 0), (1, 1)], [1, 2]), #src
+        FastDifferentiation.Node.([1.0, 2.0]), #src
+    ), #src
 ) == 3 #src
diff --git a/docs/src/2-differentiate-enzyme-model.jl b/docs/src/2-differentiate-enzyme-model.jl
index d53e2ea..c2ed938 100644
--- a/docs/src/2-differentiate-enzyme-model.jl
+++ b/docs/src/2-differentiate-enzyme-model.jl
@@ -16,13 +16,13 @@
 
 using DifferentiableMetabolism
 using AbstractFBCModels
-using Symbolics
 using ConstraintTrees
 using COBREXA
 using Tulip
 using JSONFBCModels
 import Downloads: download
 using CairoMakie
+using FastDifferentiation
 
 !isfile("e_coli_core.json") &&
     download("http://bigg.ucsd.edu/static/models/e_coli_core.json", "e_coli_core.json")
@@ -38,12 +38,10 @@ model.reactions[glc_idx]["lower_bound"] = -1000.0
 pfl_idx = first(indexin(["PFL"], rids))
 model.reactions[pfl_idx]["upper_bound"] = 0.0
 
-kcats = vcat([Symbolics.@variables $x for x in Symbol.(keys(ecoli_core_reaction_kcats))]...)
+rid_kcat = Dict(k => FastDifferentiation.Node(Symbol(k)) for (k,_) in ecoli_core_reaction_kcats)
+kcats = Symbol.(keys(ecoli_core_reaction_kcats))
 
-rid_kcat = Dict(zip(keys(ecoli_core_reaction_kcats), kcats))
-
-parameter_values =
-    Dict(kid => ecoli_core_reaction_kcats[rid] * 3.6 for (rid, kid) in rid_kcat) # change unit to k/h
+parameter_values = Dict{Symbol, Float64}()
 
 reaction_isozymes = Dict{String,Dict{String,ParameterIsozyme}}() # a mapping from reaction IDs to isozyme IDs to isozyme structs.
 float_reaction_isozymes = Dict{String,Dict{String,COBREXA.Isozyme}}() #src
@@ -52,25 +50,29 @@ for rid in AbstractFBCModels.reactions(model)
     isnothing(grrs) && continue # skip if no grr available
     haskey(ecoli_core_reaction_kcats, rid) || continue # skip if no kcat data available
     for (i, grr) in enumerate(grrs)
+
+        kcat = ecoli_core_reaction_kcats[rid] * 3.6 # change unit to k/h
+        parameter_values[Symbol(rid)] = kcat
+
         d = get!(reaction_isozymes, rid, Dict{String,ParameterIsozyme}())
-        d2 = get!(float_reaction_isozymes, rid, Dict{String,COBREXA.Isozyme}()) #src
-        d["isozyme_"*string(i)] = ParameterIsozyme(
+        d["isozyme_$i"] = ParameterIsozyme(
             gene_product_stoichiometry = Dict(grr .=> fill(1.0, size(grr))), # assume subunit stoichiometry of 1 for all isozymes
             kcat_forward = rid_kcat[rid],
             kcat_reverse = rid_kcat[rid],
         )
-        d2["isozyme_"*string(i)] = COBREXA.Isozyme( #src
+        d2 = get!(float_reaction_isozymes, rid, Dict{String,COBREXA.Isozyme}()) #src
+        d2["isozyme_$i"] = COBREXA.Isozyme( #src
             gene_product_stoichiometry = Dict(grr .=> fill(1.0, size(grr))), #src
-            kcat_forward = ecoli_core_reaction_kcats[rid] * 3.6, #src
-            kcat_reverse = ecoli_core_reaction_kcats[rid] * 3.6, #src
+            kcat_forward = kcat, #src
+            kcat_reverse = kcat, #src
         ) #src
     end
 end
 
 gene_product_molar_masses = Dict(k => v for (k, v) in ecoli_core_gene_product_masses)
 
-Symbolics.@variables capacitylimitation
-parameter_values[capacitylimitation] = 50.0 # mg enzyme/gDW
+@variables capacitylimitation
+parameter_values[:capacitylimitation] = 50.0 # mg enzyme/gDW
 
 km = build_kinetic_model(
     model;
@@ -100,15 +102,15 @@ ec_solution_cobrexa = enzyme_constrained_flux_balance_analysis( #src
 @test isapprox(ec_solution.objective, ec_solution_cobrexa.objective; atol = TEST_TOLERANCE) #src
 
 # This solution contains many inactive reactions
-sort(abs.(collect(values(ec_solution.fluxes))))
+sort(collect(ec_solution.fluxes), by=ComposedFunction(abs, last))
 
-@test any(abs.(collect(values(ec_solution.fluxes))) .≈ 0) #src
+@test any(values(ec_solution.fluxes) .≈ 0) #src
 
 # And also many inactive gene products. 
 
-sort(abs.(collect(values(ec_solution.gene_product_amounts))))
+sort(collect(ec_solution.gene_product_amounts), by=last)
 
-@test any(round.(abs.(collect(values(ec_solution.gene_product_amounts))); digits = 8) .≈ 0) #src
+@test any(isapprox.(values(ec_solution.gene_product_amounts), 0, atol=1e-8)) #src
 
 # With theory, you can show that this introduces flux variability into the
 # solution, making it non-unique, and consequently non-differentiable. To fix
@@ -121,7 +123,7 @@ flux_zero_tol = 1e-6 # these bounds make a real difference!
 gene_zero_tol = 1e-6
 
 pruned_reaction_isozymes =
-    prune_reaction_isozymes(reaction_isozymes, ec_solution, gene_zero_tol)
+    prune_reaction_isozymes(reaction_isozymes, ec_solution, flux_zero_tol)
 
 pruned_gene_product_molar_masses =
     prune_gene_product_molar_masses(gene_product_molar_masses, ec_solution, gene_zero_tol)
@@ -147,7 +149,7 @@ ec_solution2, x_vals, eq_dual_vals, ineq_dual_vals = optimized_constraints_with_
 ec_solution2
 
 # no zero fluxes
-sort(abs.(collect(values(ec_solution2.fluxes))))
+sort(collect(ec_solution.fluxes), by=ComposedFunction(abs, last))
 
 # no zero genes
 sort(abs.(collect(values(ec_solution2.gene_product_amounts))))
@@ -159,25 +161,21 @@ sort(abs.(collect(values(ec_solution2.gene_product_amounts))))
     k in intersect(keys(ec_solution.fluxes), keys(ec_solution2.fluxes)) #src
 ) #src
 
-# prune the kcats, leaving only those that are actually used
 
-pruned_kcats = [
-    begin
-        x = first(values(v))
-        isnothing(x.kcat_forward) ? x.kcat_reverse : x.kcat_forward
-    end for v in values(pruned_reaction_isozymes)
-]
-
-parameters = [capacitylimitation; kcats]
+parameters = [:capacitylimitation; kcats]
 
-sens, vids = differentiate(
+pkm_kkt, vids = differentiate_prepare_kkt(
     pkm,
     pkm.objective.value,
+    parameters
+)
+
+sens = differentiate_solution(
+    pkm_kkt,
     x_vals,
     eq_dual_vals,
     ineq_dual_vals,
     parameter_values,
-    parameters;
     scale = true, # unitless sensitivities
 )
 
diff --git a/docs/src/3-parameter-estimation.jl b/docs/src/3-parameter-estimation.jl
index ad3136c..63727e0 100644
--- a/docs/src/3-parameter-estimation.jl
+++ b/docs/src/3-parameter-estimation.jl
@@ -16,7 +16,7 @@
 
 using DifferentiableMetabolism
 using AbstractFBCModels
-using Symbolics
+using FastDifferentiation
 using ConstraintTrees
 using COBREXA
 using Clarabel
@@ -39,7 +39,7 @@ model.reactions["r2"].lower_bound = -1000.0
 
 # ![simple_model](./assets/simple_model_pruned.svg)
 
-kcats = Symbolics.@variables r3 r4
+@variables r3 r4
 
 reaction_isozymes = Dict(
     "r3" => Dict(
@@ -60,9 +60,9 @@ reaction_isozymes = Dict(
 
 gene_product_molar_masses = Dict("g1" => 20.0, "g2" => 10.0)
 
-Symbolics.@variables capacitylimitation
+@variables capacitylimitation
 
-true_parameter_values = Dict(capacitylimitation => 50.0, r3 => 2.0, r4 => 3.0)
+true_parameter_values = Dict(:capacitylimitation => 50.0, :r3 => 2.0, :r4 => 3.0)
 
 km = build_kinetic_model(
     model;
@@ -103,12 +103,17 @@ km *=
         bound = nothing,
     )
 
-estimated_parameters = Dict(capacitylimitation => 50.0, r3 => 5.0, r4 => 1.0) # initial values
-parameters = [r3, r4] # will differentiate against these two parameters
+estimated_parameters = Dict(:capacitylimitation => 50.0, :r3 => 5.0, :r4 => 1.0) # initial values
 η = 0.1 # learning rate
 
 losses = Float64[]
 
+kmKKT, vids = differentiate_prepare_kkt(
+    km,
+    km.loss.value,
+    [:r3, :r4, :capacitylimitation],
+)
+
 for k = 1:150
 
     sol2, x_vals, eq_dual_vals, ineq_dual_vals = optimized_constraints_with_parameters(
@@ -122,16 +127,14 @@ for k = 1:150
 
     push!(losses, sol2.loss)
 
-    sens, vids = differentiate(
-        km,
-        km.loss.value,
+    sens = differentiate_solution(
+        kmKKT,
         x_vals,
         eq_dual_vals,
         ineq_dual_vals,
         estimated_parameters,
-        parameters;
     )
-    measured_idxs = [1, 3, 11, 12]
+    measured_idxs = [1, 3, 12, 11]
 
     x = [
         sol2.fluxes.r1,
@@ -143,8 +146,8 @@ for k = 1:150
     dL_dx = x - measured # derivative of loss function with respect to optimization variables
     dL_dkcats = sens[measured_idxs, :]' * dL_dx # derivative of loss function with respect to parameters
 
-    estimated_parameters[r3] -= η * dL_dkcats[1]
-    estimated_parameters[r4] -= η * dL_dkcats[2]
+    estimated_parameters[:r3] -= η * dL_dkcats[1]
+    estimated_parameters[:r4] -= η * dL_dkcats[2]
 end
 
 lines(losses; axis = (xlabel = "Iterations", ylabel = "L2 loss"))
@@ -156,7 +159,6 @@ estimated_parameters
 
 true_parameter_values
 
-@test abs(estimated_parameters[r3] - true_parameter_values[r3]) <= 0.1 #src
-@test abs(estimated_parameters[r4] - true_parameter_values[r4]) <= 0.1 #src
-@test abs(estimated_parameters[r4] - true_parameter_values[r4]) <= 0.1 #src
+@test abs(estimated_parameters[:r3] - true_parameter_values[:r3]) <= 0.1 #src
+@test abs(estimated_parameters[:r4] - true_parameter_values[:r4]) <= 0.1 #src
 @test all(losses[2:end] .<= losses[1:end-1]) #src
diff --git a/docs/src/reference.md b/docs/src/reference.md
index 2aa9f73..40a4ca0 100644
--- a/docs/src/reference.md
+++ b/docs/src/reference.md
@@ -76,9 +76,9 @@ Modules = [DifferentiableMetabolism]
 Pages = ["src/get_constraints.jl"]
 ```
 
-### Symbolics extensions
+### Expression evaluation utilities
 
 ```@autodocs
 Modules = [DifferentiableMetabolism]
-Pages = ["src/symbolics.jl"]
+Pages = ["src/substitute.jl"]
 ```
diff --git a/src/.DS_Store b/src/.DS_Store
deleted file mode 100644
index 1cbfc23..0000000
Binary files a/src/.DS_Store and /dev/null differ
diff --git a/src/differentiate.jl b/src/differentiate.jl
index 065871d..4d9c6a6 100644
--- a/src/differentiate.jl
+++ b/src/differentiate.jl
@@ -51,6 +51,8 @@ function findall_indeps_qr(A)
     t.pcol[1:max_lin_indep_columns] # undo permumation
 end
 
+export differentiate_prepare_kkt
+
 """
 $(TYPEDSIGNATURES)
 
@@ -124,6 +126,20 @@ function differentiate_prepare_kkt(
     A = FastDifferentiation.sparse_jacobian(kkt_eqns, [xs; eq_duals; ineq_duals])
     B = FastDifferentiation.sparse_jacobian(kkt_eqns, FastDifferentiation.Node.(parameters))
 
+    return (A, B, xs, eq_duals, ineq_duals, parameters), variable_order(m)
+end
+
+export differentiate_prepare_kkt
+
+function differentiate_solution(
+    (A, B, xs, eq_duals, ineq_duals, parameters),
+    x_vals::Vector{Float64},
+    eq_dual_vals::Vector{Float64},
+    ineq_dual_vals::Vector{Float64},
+    parameter_values::Dict{Symbol,Float64};
+    scale = false, # scale sensitivities
+)
+
     # symbolic values at the optimal solution incl parameters
     syms_to_vals = merge(
         Dict(zip((x.node_value for x in [xs; eq_duals; ineq_duals]), [x_vals; eq_dual_vals; ineq_dual_vals])),
@@ -152,9 +168,9 @@ function differentiate_prepare_kkt(
 
     # get primal variable sensitivities only
     if scale
-        ([parameter_values[p] for p=parameters]' .* c ./ x_vals)
+        ([parameter_values[p] for p in parameters]' .* c[1:length(xs), :] ./ x_vals)
     else
-        c
+        c[1:length(xs),:]
     end
 end
 
diff --git a/src/parameter_linearvalue.jl b/src/parameter_linearvalue.jl
index 8c8b2a5..06cedca 100644
--- a/src/parameter_linearvalue.jl
+++ b/src/parameter_linearvalue.jl
@@ -37,7 +37,7 @@ $(TYPEDFIELDS)
     weights::Vector{Expression}
 end
 
-# convert all weights to Nums
+# convert all weights to Expressions
 ParameterLinearValue(idxs::Vector{Int}, weights::Vector{Union{Int,Float64}}) =
     ParameterLinearValue(idxs, convert.(Expression, weights))
 
diff --git a/src/parameter_promotion.jl b/src/parameter_promotion.jl
index 95eff05..352939e 100644
--- a/src/parameter_promotion.jl
+++ b/src/parameter_promotion.jl
@@ -18,7 +18,7 @@ limitations under the License.
 
 # Promote LinearValue to ParameterLinearValue if it interacts with a parameter
 
-# Nums and LinearValues
+# Expressions and LinearValues
 Base.:+(a::Expression, b::ConstraintTrees.LinearValue) = b + a
 Base.:+(a::ConstraintTrees.LinearValue, b::Expression) = a + ParameterLinearValue(b)
 
@@ -43,7 +43,7 @@ Base.:-(a::ParameterLinearValue, b::ConstraintTrees.LinearValue) =
 
 # Promote QuadraticValue to ParameterQuadraticValue if it interacts with a parameter
 
-# Nums and QuadraticValues
+# Expressions and QuadraticValues
 Base.:+(a::Expression, b::ConstraintTrees.QuadraticValue) = b + a
 Base.:+(a::ConstraintTrees.QuadraticValue, b::Expression) =
     a + ParameterQuadraticValue(b)