diff --git a/docs/src/examples/02a-optimizer-parameters.jl b/docs/src/examples/02a-optimizer-parameters.jl
index 1a81f20cc..0bc355c99 100644
--- a/docs/src/examples/02a-optimizer-parameters.jl
+++ b/docs/src/examples/02a-optimizer-parameters.jl
@@ -17,7 +17,7 @@
 # # Changing optimizer parameters
 #
 # Many optimizers require fine-tuning to produce best results. You can pass in
-# additional optimizer settings via the `modifications` parameter of
+# additional optimizer settings via the `settings` parameter of
 # [`flux_balance_analysis`](@ref). These include e.g.
 #
 # - [`set_optimizer_attribute`](@ref) (typically allowing you to tune e.g.
@@ -47,7 +47,7 @@ model = load_model("e_coli_core.json")
 solution = flux_balance_analysis(
     model,
     Tulip.Optimizer;
-    modifications = [silence, set_optimizer_attribute("IPM_IterationsLimit", 1000)],
+    settings = [silence, set_optimizer_attribute("IPM_IterationsLimit", 1000)],
 )
 
 @test !isnothing(solution) #src
@@ -60,7 +60,7 @@ solution = flux_balance_analysis(
 solution = flux_balance_analysis(
     model,
     Tulip.Optimizer;
-    modifications = [set_optimizer_attribute("IPM_IterationsLimit", 2)],
+    settings = [set_optimizer_attribute("IPM_IterationsLimit", 2)],
 )
 
 println(solution)
diff --git a/docs/src/examples/02c-constraint-modifications.jl b/docs/src/examples/02c-constraint-modifications.jl
index 07f0f690e..9afaa534f 100644
--- a/docs/src/examples/02c-constraint-modifications.jl
+++ b/docs/src/examples/02c-constraint-modifications.jl
@@ -62,7 +62,7 @@ vt = optimized_constraints(
     forced_mixed_fermentation,
     objective = forced_mixed_fermentation.objective.value,
     optimizer = Tulip.Optimizer,
-    modifications = [silence],
+    settings = [silence],
 )
 
 @test isapprox(vt.objective, 0.6337, atol = TEST_TOLERANCE) #src
@@ -78,7 +78,7 @@ vt = optimized_constraints(
     ctmodel,
     objective = ctmodel.objective.value,
     optimizer = Tulip.Optimizer,
-    modifications = [silence],
+    settings = [silence],
 )
 
 @test isnothing(vt) #src
@@ -90,7 +90,7 @@ vt = optimized_constraints(
     ctmodel,
     objective = ctmodel.objective.value,
     optimizer = Tulip.Optimizer,
-    modifications = [silence],
+    settings = [silence],
 )
 
 @test isapprox(vt.objective, 0.8739, atol = TEST_TOLERANCE) #src
diff --git a/docs/src/examples/03-parsimonious-flux-balance.jl b/docs/src/examples/03-parsimonious-flux-balance.jl
index d69ad5a07..eee4439ec 100644
--- a/docs/src/examples/03-parsimonious-flux-balance.jl
+++ b/docs/src/examples/03-parsimonious-flux-balance.jl
@@ -40,11 +40,11 @@ model = load_model("e_coli_core.json") # load the model
 # Use the convenience function to run standard pFBA on
 
 vt =
-    parsimonious_flux_balance_analysis(model, Clarabel.Optimizer; modifications = [silence])
+    parsimonious_flux_balance_analysis(model, Clarabel.Optimizer; settings = [silence])
 
 # Or use the piping functionality
 
-model |> parsimonious_flux_balance_analysis(Clarabel.Optimizer; modifications = [silence])
+model |> parsimonious_flux_balance_analysis(Clarabel.Optimizer; settings = [silence])
 
 @test isapprox(vt.objective, 0.87392; atol = TEST_TOLERANCE) #src
 @test sum(x^2 for x in values(vt.fluxes)) < 15000 #src
@@ -82,7 +82,7 @@ vt = minimize_metabolic_adjustment(model, ref_sol, Gurobi.Optimizer)
 # Or use the piping functionality
 
 model |>
-minimize_metabolic_adjustment(ref_sol, Clarabel.Optimizer; modifications = [silence])
+minimize_metabolic_adjustment(ref_sol, Clarabel.Optimizer; settings = [silence])
 
 @test isapprox(vt.:momaobjective, 0.81580806; atol = TEST_TOLERANCE) #src
 
diff --git a/docs/src/examples/05-enzyme-constrained-models.jl b/docs/src/examples/05-enzyme-constrained-models.jl
index a65c178e9..caad58073 100644
--- a/docs/src/examples/05-enzyme-constrained-models.jl
+++ b/docs/src/examples/05-enzyme-constrained-models.jl
@@ -120,7 +120,7 @@ ec_solution = enzyme_constrained_flux_balance_analysis(
     reaction_isozymes,
     gene_molar_masses,
     [("total_proteome_bound", A.genes(model), total_enzyme_capacity)];
-    modifications = [set_optimizer_attribute("IPM_IterationsLimit", 10_000)],
+    settings = [set_optimizer_attribute("IPM_IterationsLimit", 10_000)],
     unconstrain_reactions = ["EX_glc__D_e"],
     optimizer = Tulip.Optimizer,
 )
@@ -173,5 +173,5 @@ ec_solution = optimized_constraints(
     m;
     objective = m.objective.value,
     optimizer = Tulip.Optimizer,
-    modifications = [set_optimizer_attribute("IPM_IterationsLimit", 10_000)],
+    settings = [set_optimizer_attribute("IPM_IterationsLimit", 10_000)],
 )
diff --git a/docs/src/examples/06-thermodynamic-models.jl b/docs/src/examples/06-thermodynamic-models.jl
index 603192272..369492bc5 100644
--- a/docs/src/examples/06-thermodynamic-models.jl
+++ b/docs/src/examples/06-thermodynamic-models.jl
@@ -116,7 +116,7 @@ mmdf_solution = max_min_driving_force_analysis(
     concentration_ub = 1e-1, # M
     T = 298.15, # Kelvin
     R = 8.31446261815324e-3, # kJ/K/mol
-    modifications = [set_optimizer_attribute("IPM_IterationsLimit", 1_000)],
+    settings = [set_optimizer_attribute("IPM_IterationsLimit", 1_000)],
     optimizer = Tulip.Optimizer,
 )
 
@@ -150,5 +150,5 @@ mmdf_solution = optimized_constraints(
     m;
     objective = m.max_min_driving_force.value,
     optimizer = Tulip.Optimizer,
-    modifications = [set_optimizer_attribute("IPM_IterationsLimit", 1_000)],
+    settings = [set_optimizer_attribute("IPM_IterationsLimit", 1_000)],
 )
diff --git a/docs/src/examples/08-community-models.jl b/docs/src/examples/08-community-models.jl
index 55a1fc389..a044b7ecc 100644
--- a/docs/src/examples/08-community-models.jl
+++ b/docs/src/examples/08-community-models.jl
@@ -110,7 +110,7 @@ sol = optimized_constraints(
     m;
     objective = m.objective.value,
     optimizer = Tulip.Optimizer,
-    modifications = [set_optimizer_attribute("IPM_IterationsLimit", 1000)],
+    settings = [set_optimizer_attribute("IPM_IterationsLimit", 1000)],
 )
 
 @test isapprox(sol.:objective, 0.66686196344, atol = TEST_TOLERANCE) #src
@@ -127,7 +127,7 @@ sol = optimized_constraints(
     m;
     objective = m.objective.value,
     optimizer = Tulip.Optimizer,
-    modifications = [set_optimizer_attribute("IPM_IterationsLimit", 1000)],
+    settings = [set_optimizer_attribute("IPM_IterationsLimit", 1000)],
 )
 
 
diff --git a/src/COBREXA.jl b/src/COBREXA.jl
index 820186124..43748a909 100644
--- a/src/COBREXA.jl
+++ b/src/COBREXA.jl
@@ -69,7 +69,7 @@ include("frontend/enzyme_constrained.jl")
 
 # utilities
 include("misc/bounds.jl")
-include("misc/modifications.jl")
+include("misc/settings.jl")
 include("misc/utils.jl")
 
 end # module COBREXA
diff --git a/src/frontend/balance.jl b/src/frontend/balance.jl
index 82e28e3ae..624f584eb 100644
--- a/src/frontend/balance.jl
+++ b/src/frontend/balance.jl
@@ -39,7 +39,7 @@ $(TYPEDSIGNATURES)
 
 Pipe-able overload of [`flux_balance_analysis`](@ref).
 """
-flux_balance_analysis(optimizer; modifications = []) =
-    m -> flux_balance_analysis(m, optimizer; modifications)
+flux_balance_analysis(optimizer; settings = []) =
+    m -> flux_balance_analysis(m, optimizer; settings)
 
 export flux_balance_analysis
diff --git a/src/frontend/enzyme_constrained.jl b/src/frontend/enzyme_constrained.jl
index 16b7a0a1b..ab0779c02 100644
--- a/src/frontend/enzyme_constrained.jl
+++ b/src/frontend/enzyme_constrained.jl
@@ -25,7 +25,7 @@ The latter is a vector of tuples, where each tuple represents a distinct bound
 as `(bound_id, genes_in_bound, protein_mass_bound)`. Typically, `model` has
 bounded exchange reactions, which are unnecessary in enzyme constrained models.
 Unbound these reactions by listing their IDs in `unconstrain_reactions`, which
-makes them reversible. Optimization `modifications` are directly forwarded.
+makes them reversible. Optimization `settings` are directly forwarded.
 
 In the event that your model requires more complex build steps, consider
 constructing it manually by using [`add_enzyme_constraints!`](@ref).
@@ -37,7 +37,7 @@ function enzyme_constrained_flux_balance_analysis(
     capacity_limitations::Vector{Tuple{String,Vector{String},Float64}};
     optimizer,
     unconstrain_reactions = String[],
-    modifications = [],
+    settings = [],
 )
     m = fbc_model_constraints(model)
 
@@ -56,7 +56,7 @@ function enzyme_constrained_flux_balance_analysis(
         m.fluxes[rid].bound = C.Between(-1000.0, 1000.0)
     end
 
-    optimized_constraints(m; objective = m.objective.value, optimizer, modifications)
+    optimized_constraints(m; objective = m.objective.value, optimizer, settings)
 end
 
 export enzyme_constrained_flux_balance_analysis
diff --git a/src/frontend/loopless.jl b/src/frontend/loopless.jl
index 6d1c49ace..0c1a8bc0c 100644
--- a/src/frontend/loopless.jl
+++ b/src/frontend/loopless.jl
@@ -41,7 +41,7 @@ function loopless_flux_balance_analysis(
     model;
     max_flux_bound = 1000.0, # needs to be an order of magnitude bigger, big M method heuristic
     strict_inequality_tolerance = 1.0, # heuristic from paper
-    modifications = [],
+    settings = [],
     optimizer,
 )
 
@@ -68,7 +68,7 @@ function loopless_flux_balance_analysis(
     )
 
     # solve
-    optimized_constraints(m; objective = m.objective.value, optimizer, modifications)
+    optimized_constraints(m; objective = m.objective.value, optimizer, settings)
 end
 
 export loopless_flux_balance_analysis
diff --git a/src/frontend/mmdf.jl b/src/frontend/mmdf.jl
index 1df1251f4..55be1fe01 100644
--- a/src/frontend/mmdf.jl
+++ b/src/frontend/mmdf.jl
@@ -67,7 +67,7 @@ supplied).
 `T` and `R` can be specified in the corresponding units; defaults are K and
 kJ/K/mol. The unit of metabolite concentrations is typically molar, and the ΔG⁰s
 have units of kJ/mol. Other units can be used, as long as they are consistent.
-As usual, optimizer settings can be changed with `modifications`.
+As usual, optimizer settings can be changed with `settings`.
 """
 function max_min_driving_force_analysis(
     model::A.AbstractFBCModel,
@@ -82,7 +82,7 @@ function max_min_driving_force_analysis(
     T = 298.15, # Kelvin
     R = 8.31446261815324e-3, # kJ/K/mol
     ignore_reaction_ids = String[],
-    modifications = [],
+    settings = [],
     optimizer,
 )
     m = build_max_min_driving_force_model(
@@ -113,7 +113,7 @@ function max_min_driving_force_analysis(
         m;
         objective = m.max_min_driving_force.value,
         optimizer,
-        modifications,
+        settings,
     )
 end
 
diff --git a/src/frontend/moma.jl b/src/frontend/moma.jl
index b1af4944e..b3acc69f0 100644
--- a/src/frontend/moma.jl
+++ b/src/frontend/moma.jl
@@ -57,7 +57,7 @@ A slightly easier-to-use version of [`minimization_of_metabolic_adjustment_analy
 computes the reference flux as the optimal solution of the
 [`reference_model`](@ref). The reference flux is calculated using
 `reference_optimizer` and `reference_modifications`, which default to the
-`optimizer` and `modifications`.
+`optimizer` and `settings`.
 
 Leftover arguments are passed to the overload of
 [`minimization_of_metabolic_adjustment_analysis`](@ref) that accepts the reference flux
@@ -68,15 +68,15 @@ function minimization_of_metabolic_adjustment_analysis(
     reference_model::A.AbstractFBCModel,
     optimizer;
     reference_optimizer = optimizer,
-    modifications = [],
-    reference_modifications = modifications,
+    settings = [],
+    reference_settings = settings,
     kwargs...,
 )
     reference_constraints = fbc_model_constraints(reference_model)
     reference_fluxes = optimized_constraints(
         reference_constraints;
         optimizer = reference_optimizer,
-        modifications = reference_modifications,
+        settings = reference_settings,
         output = reference_constraints.fluxes,
     )
     isnothing(reference_fluxes) && return nothing
@@ -84,7 +84,7 @@ function minimization_of_metabolic_adjustment_analysis(
         model,
         reference_fluxes,
         optimizer;
-        modifications,
+        settings,
         kwargs...,
     )
 end
diff --git a/src/frontend/parsimonious.jl b/src/frontend/parsimonious.jl
index c76543082..1f7d732f5 100644
--- a/src/frontend/parsimonious.jl
+++ b/src/frontend/parsimonious.jl
@@ -20,7 +20,7 @@ $(TYPEDSIGNATURES)
 Optimize the system of `constraints` to get the optimal `objective` value. Then
 try to find a "parsimonious" solution with the same `objective` value, which
 optimizes the `parsimonious_objective` (possibly also switching optimization
-sense, optimizer, and adding more modifications).
+sense, optimizer, and adding more settings).
 
 For efficiency, everything is performed on a single instance of JuMP model.
 
@@ -30,7 +30,7 @@ in [`parsimonious_flux_balance`](@ref).
 function parsimonious_optimized_constraints(
     constraints::C.ConstraintTreeElem;
     objective::C.Value,
-    modifications = [],
+    settings = [],
     parsimonious_objective::C.Value,
     parsimonious_optimizer = nothing,
     parsimonious_sense = J.MIN_SENSE,
@@ -42,7 +42,7 @@ function parsimonious_optimized_constraints(
 
     # first solve the optimization problem with the original objective
     om = optimization_model(constraints; objective, kwargs...)
-    for m in modifications
+    for m in settings
         m(om)
     end
     J.optimize!(om)
diff --git a/src/misc/modifications.jl b/src/misc/settings.jl
similarity index 100%
rename from src/misc/modifications.jl
rename to src/misc/settings.jl
diff --git a/src/solver.jl b/src/solver.jl
index d9f4f85c3..041c677cf 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -104,7 +104,7 @@ export Feasible
 $(TYPEDSIGNATURES)
 
 Make an JuMP model out of `constraints` using [`optimization_model`](@ref)
-(most arguments are forwarded there), then apply the `modifications`, optimize
+(most arguments are forwarded there), then apply the `settings`, optimize
 the model, and return either `nothing` if the optimization failed, or `output`
 substituted with the solved values (`output` defaults to `constraints`.
 
@@ -113,12 +113,12 @@ For a "nice" version for simpler finding of metabolic model optima, use
 """
 function optimized_constraints(
     constraints::C.ConstraintTreeElem;
-    modifications = [],
+    settings = [],
     output = constraints,
     kwargs...,
 )
     om = optimization_model(constraints; kwargs...)
-    for m in modifications
+    for m in settings
         m(om)
     end
     J.optimize!(om)