diff --git a/src/analysis/fba.jl b/src/analysis/fba.jl
index 2e15bf838e2017034e8883ab35aa34d1df4d427a..363a868411976395d6b7ec840e12619b122264f3 100644
--- a/src/analysis/fba.jl
+++ b/src/analysis/fba.jl
@@ -71,8 +71,9 @@ function fba(
     sense = MOI.MAX_SENSE,
 )
     # get core optimization problem
-    cbm, v, mb, lbcons, ubcons = make_optimization_model(model, optimizer, sense = sense)
-
+    cbm = make_optimization_model(model, optimizer, sense = sense)
+    v = cbm[:x] # fluxes
+    
     # modify core optimization problem according to user specifications
     if !isempty(solver_attributes) # set other attributes
         for (k, val) in solver_attributes
@@ -83,7 +84,7 @@ function fba(
     # set additional constraints
     for (rxnid, con) in constraints
         ind = model.reactions[findfirst(model.reactions, rxnid)]
-        set_bound(ind, lbcons, ubcons; lb = con[1], ub = con[2])
+        set_bound(ind, cbm; lb = con[1], ub = con[2])
     end
 
     # if an objective function is supplied, modify the default objective
diff --git a/src/analysis/fva.jl b/src/analysis/fva.jl
index 0b38dc87179e3d82db28c96b2f5abcd4db8a0f6c..e1c5a575c1470b3341423f9d3095d1cd532efda9 100644
--- a/src/analysis/fva.jl
+++ b/src/analysis/fva.jl
@@ -37,14 +37,15 @@ function flux_variability_analysis(
         throw(DomainError(reactions, "Index exceeds number of reactions."))
     end
 
-    (optimization_model, x0) = flux_balance_analysis(model, optimizer)
+    optimization_model = flux_balance_analysis(model, optimizer)
     Z0 = COBREXA.JuMP.objective_value(optimization_model)
     optimization_model = nothing # we won't need this one anymore, so free the memory
 
     # store a JuMP optimization model at all workers
     save_model = :(
         begin
-            optmodel, x = COBREXA.make_optimization_model($model, $optimizer)
+            optmodel = COBREXA.make_optimization_model($model, $optimizer)
+            x = optmodel[:x]
             COBREXA._FVA_add_constraint(optmodel, $(objective(model)), x, $Z0, $gamma)
             optmodel
         end
@@ -140,7 +141,8 @@ function fva(
     constraints = Dict{String,Tuple{Float64,Float64}}(),
     sense = MOI.MAX_SENSE,
 )
-    cbm, v, mb, lbcons, ubcons = make_optimization_model(model, optimizer, sense = sense)
+    cbm = make_optimization_model(model, optimizer, sense = sense)
+    v = cbm[:x]
 
     if !isempty(solver_attributes) # set other attributes
         for (k, v) in solver_attributes
@@ -151,7 +153,7 @@ function fva(
     # set additional constraints
     for (rxnid, con) in constraints
         ind = model.reactions[findfirst(model.reactions, rxnid)]
-        set_bound(ind, lbcons, ubcons; lb = con[1], ub = con[2])
+        set_bound(ind, cbm; lb = con[1], ub = con[2])
     end
 
     # if an objective function is supplied, modify the default objective
diff --git a/src/analysis/pfba.jl b/src/analysis/pfba.jl
index 1e2afc4139a4494610cda97493d8da0643a5a941..8402d3f4da2052a72ddc9f9003a211a1d44c4047 100644
--- a/src/analysis/pfba.jl
+++ b/src/analysis/pfba.jl
@@ -33,16 +33,16 @@ function pfba(
 
 
     if typeof(optimizer) <: AbstractArray # choose optimizer
-        cbm, v, mb, lbcons, ubcons =
-            make_optimization_model(model, optimizer[1], sense = sense)
+        cbm = make_optimization_model(model, optimizer[1], sense = sense)
+        v = cbm[:x]
         if !isempty(solver_attributes["opt1"]) # set other attributes
             for (k, v) in solver_attributes["opt1"]
                 set_optimizer_attribute(cbm, k, v)
             end
         end
     else # singe optimizer
-        cbm, v, mb, lbcons, ubcons =
-            make_optimization_model(model, optimizer, sense = sense)
+        cbm = make_optimization_model(model, optimizer, sense = sense)
+        v = cbm[:x]
         if !isempty(solver_attributes) # set other attributes
             for (k, v) in solver_attributes
                 set_optimizer_attribute(cbm, k, v)
@@ -53,7 +53,7 @@ function pfba(
     # set additional constraints
     for (rxnid, con) in constraints
         ind = model.reactions[findfirst(model.reactions, rxnid)]
-        set_bound(ind, lbcons, ubcons; lb = con[1], ub = con[2])
+        set_bound(ind, cbm; lb = con[1], ub = con[2])
     end
 
     # check if default objective should be used
diff --git a/src/base/solver.jl b/src/base/solver.jl
index 178654463f73c57a05844438ca73c0cdf8674620..b8d16887170516b01d66a9bc8d4e31c08e5f8a66 100644
--- a/src/base/solver.jl
+++ b/src/base/solver.jl
@@ -13,11 +13,11 @@ function make_optimization_model(
     optimization_model = COBREXA.JuMP.Model(optimizer)
     @variable(optimization_model, x[i = 1:n])
     @objective(optimization_model, sense, objective(model)' * x)
-    mb = @constraint(optimization_model, mb, stoichiometry(model) * x .== balance(model)) # mass balance
-    lbs = @constraint(optimization_model, lbs, xl .<= x) # lower bounds
-    ubs = @constraint(optimization_model, ubs, x .<= xu) # upper bounds
+    @constraint(optimization_model, mb, stoichiometry(model) * x .== balance(model)) # mass balance
+    @constraint(optimization_model, lbs, xl .<= x) # lower bounds
+    @constraint(optimization_model, ubs, x .<= xu) # upper bounds
 
-    return optimization_model, x, mb, lbs, ubs
+    return optimization_model
 end
 
 """
@@ -28,8 +28,7 @@ function optimize_model(
     optimizer;
     sense = MOI.MIN_SENSE,
 ) where {LM<:MetabolicModel}
-    optimization_model, x, _, _, _ =
-        make_optimization_model(model, optimizer; sense = sense)
-        COBREXA.JuMP.optimize!(optimization_model)
-    return (optimization_model, x)
+    optimization_model = make_optimization_model(model, optimizer; sense = sense)
+    COBREXA.JuMP.optimize!(optimization_model)
+    return optimization_model
 end
diff --git a/src/base/utilities.jl b/src/base/utilities.jl
index 6ea40a1767aae5e13e9b23f6a48d0c07c804d242..ac473dd5677cb0e740583539a7939ae4fabd251d 100644
--- a/src/base/utilities.jl
+++ b/src/base/utilities.jl
@@ -21,16 +21,16 @@ Base.copy(model::CoupledLinearModel) =
 
 
 """
-    set_bound(index, ubconstaintref, lbconstaintref; ub=1000, lb=-1000)
+    set_bound(index, optimization_model; ub=1000, lb=-1000)
 Helper function to set the bounds of variables.
 The JuMP `set_normalized_rhs` function is a little confusing, so this function simplifies setting
-constraints. Just supply the constraint `index` and the bound objects (`ubconstaintref`, `lbconstaintref`) and they will be set to `ub` and `lb`.
+constraints. Just supply the constraint `index` and the model and they will be set to `ub` and `lb`.
 """
-function set_bound(vind, lbs, ubs; ub = 1000, lb = -1000)
+function set_bound(vind, opt_model; ub = 1000, lb = -1000)
     if lb <= 0
-        set_normalized_rhs(lbs[vind], abs(lb))
+        set_normalized_rhs(opt_model[:lbs][vind], abs(lb))
     else
-        set_normalized_rhs(lbs[vind], -abs(lb))
+        set_normalized_rhs(opt_model[:lbs][vind], -abs(lb))
     end
-    set_normalized_rhs(ubs[vind], ub)
+    set_normalized_rhs(opt_model[:ubs][vind], ub)
 end