diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6e50d57..5c19cb2 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Unreleased
 
+* `get_conc_control_coeffs` takes the control coefficient by finite difference by default, as GECKO's `getConcControlCoeffs` does: it raises one usage bound by `fold_change` (2), re-solves and divides the rise of the optimum by the rise of the bound. The reduced-cost path, which reads the dual of the usage reaction and is unreliable at a degenerate vertex, is `method="shadow_price"`. The finite-difference path needs one solve per protein, so each iteration of `flexibilize_enz_concs`, which calls it, takes minutes on a genome-scale ecModel.
 * `get_subset_ec_model` works on the scipy solver. It raised `LookupError: Constraint prot_P1 ... not in solver` at the next `optimize` when it had removed two or more constraints: optlang's scipy interface calls the generic removal inside its per-constraint loop, with the whole list each time, so the second pass finds the constraints gone. The sparse scipy model removes constraints correctly, and the function removes genes, reactions and metabolites on glpk and converts the subset back to scipy, because optlang rewrites every constraint expression for each removed variable (65 minutes for the HT29-GEM subset of the light Human-GEM ecModel). `geckopy.solver_support.edit_on_fast_solver` is the context manager that does this for any structural edit of a scipy model.
 * `switch_to_scipy` and so `make_ec_model`, `load_conventional_gem` and `load_ec_model` on scipy create each constraint with its whole expression instead of setting the coefficients one at a time, which made the scipy interface remove and add the rows of a constraint for every call (the light Human-GEM ecModel took 2 h 17 min to build on scipy against 1.5 min on Gurobi). `apply_kcat_constraints` edits the model on glpk when it updates 100 or more ec reactions of a scipy model, for the same reason.
 * `ec_fseof` follows `ecFSEOF.m`. The target flux is fixed at each level of the scan (`equality=True` of the raven-toolbox `fseof`, which `ecFSEOF.m` sets with `setParam(model, 'eq', ...)`) instead of lower-bounded. `targets` lists the reaction targets first and then the transport targets (MATLAB's `rxnTargets` and `transportTargets`), each group sorted by score, and has a new boolean column `transport`; a reaction is a transport target when a metabolite name occurs more than once among its metabolites, not counting `prot_*`. This changes the row order of `targets` (and of `amplification`, `knockout` and `gene_targets`, which are derived from it) for callers that take the top N: transport targets now come after all reaction targets.
diff --git a/src/geckopy/limit_proteins/get_conc_control_coeffs.py b/src/geckopy/limit_proteins/get_conc_control_coeffs.py
index 280634b..9b09fda 100644
--- a/src/geckopy/limit_proteins/get_conc_control_coeffs.py
+++ b/src/geckopy/limit_proteins/get_conc_control_coeffs.py
@@ -47,6 +47,7 @@ def get_conc_control_coeffs(
     proteins: Optional[list[str]] = None,
     fold_change: float = 2.0,
     limit: float = 0.0,
+    method: str = "finite_difference",
 ) -> tuple[np.ndarray, np.ndarray]:
     """Compute control coefficients for enzyme concentrations.
 
@@ -56,7 +57,14 @@ def get_conc_control_coeffs(
     ``limit`` times their upper bound are skipped (no sensitivity to
     expect there).
 
-    The default path reads, from a single LP solve, the reduced cost of
+    The default ``method="finite_difference"`` raises the upper bound of one usage reaction
+    by ``fold_change``, re-solves and divides the rise of the optimum by the rise of the
+    bound, one solve per analysed protein (several minutes on a genome-scale ecModel). A
+    rise of at most 1e-10 counts as no sensitivity. This is what GECKO's
+    ``getConcControlCoeffs`` does, and it is a re-solved change of the objective, which
+    stays meaningful when the LP is degenerate and its duals are not.
+
+    ``method="shadow_price"`` reads, from a single LP solve, the reduced cost of
     each analysed usage reaction's forward variable. That is the local
     marginal d(objective)/d(usage_ub) of the enzyme's own bound: it is
     zero for an enzyme whose bound does not bind, also when the shared
@@ -73,9 +81,9 @@ def get_conc_control_coeffs(
     read from the bound marginal of the last ``linprog`` result, since optlang does not
     expose variable duals there. When the active solver gives no usable variable duals
     (a missing or NaN value, or a scipy result without bound marginals) the function falls
-    back to the per-protein finite-difference loop, which needs one solve per protein
-    (several minutes on a genome-scale ecModel), and logs once that it has done so. The
-    fold-change probe size is only used in that fallback.
+    back to the finite-difference loop and logs once that it has done so. The two methods
+    agree where raising the bound does not cross a change of the active set and differ
+    where it does.
 
     Ported from GECKO MATLAB:
     src/geckomat/limit_proteins/getConcControlCoeffs.m.
@@ -89,12 +97,14 @@ def get_conc_control_coeffs(
     proteins
         UniProt IDs to analyse. Defaults to ``model.ec.enzymes``.
     fold_change
-        Probe size for the finite-difference fallback (default 2.0).
-        Ignored on the shadow-price path.
+        Factor by which the finite-difference path raises a usage upper bound
+        (default 2.0). Used by ``method="shadow_price"`` only when it falls back.
     limit
         Skip proteins whose usage flux is at most ``limit`` times
         their upper bound. Default 0.0 includes any protein with
         non-zero flux.
+    method
+        ``"finite_difference"`` (default) or ``"shadow_price"``.
 
     Returns
     -------
@@ -105,13 +115,18 @@ def get_conc_control_coeffs(
         Float array of length ``len(proteins)``; 0.0 where the
         protein was either skipped or showed no sensitivity.
     """
+    if method not in ("finite_difference", "shadow_price"):
+        raise ValueError(
+            f"get_conc_control_coeffs: method must be 'finite_difference' or "
+            f"'shadow_price', got {method!r}"
+        )
     if proteins is None:
         proteins = list(model.ec.enzymes)
-
     n = len(proteins)
     if n == 0:
         return np.zeros(0, dtype=bool), np.zeros(0, dtype=float)
-
+    if method == "finite_difference":
+        return _finite_difference_coeffs(model, proteins, fold_change, limit)
     if _solver_supports_duals(model):
         return _shadow_price_coeffs(model, proteins, limit, fold_change)
     _report_finite_difference_fallback(model, "does not expose LP duals")
@@ -240,10 +255,10 @@ def _finite_difference_coeffs(
     fold_change: float,
     limit: float,
 ) -> tuple[np.ndarray, np.ndarray]:
-    """Per-protein finite-difference fallback (one solve per protein).
+    """Per-protein finite-difference coefficients (one solve per protein).
 
-    Same contract as :func:`get_conc_control_coeffs`. Used when the solver
-    does not expose usable variable duals.
+    Same contract as :func:`get_conc_control_coeffs`. The default method, and the
+    fallback of the shadow-price method when the solver does not expose usable variable duals.
     """
     n = len(proteins)
     enz_mask = np.zeros(n, dtype=bool)
diff --git a/tests/test_get_conc_control_coeffs.py b/tests/test_get_conc_control_coeffs.py
index 05454bb..c555b48 100644
--- a/tests/test_get_conc_control_coeffs.py
+++ b/tests/test_get_conc_control_coeffs.py
@@ -438,13 +438,13 @@ def test_finite_diff_underestimates_when_step_crosses_breakpoint():
 # finite-difference where the solver gives no usable variable duals.
 # --------------------------------------------------------------------------- #
 
-def test_dispatcher_uses_shadow_price_on_default_solver():
-    """On glpk (cobra default), the public function takes the shadow-price
-    path. Verified by matching the shadow-price helper exactly (the
+def test_dispatcher_uses_shadow_price_when_asked():
+    """With method="shadow_price" on glpk (cobra default), the public function takes the
+    shadow-price path. Verified by matching the shadow-price helper exactly (the
     finite-difference would diverge on a breakpoint scenario)."""
     model = _build_toy_ec_model(enzyme_ub=8.0)
     model.reactions.get_by_id("EX_A").lower_bound = -10.0  # breakpoint setup
-    pm, pc = get_conc_control_coeffs(model, fold_change=2.0)
+    pm, pc = get_conc_control_coeffs(model, fold_change=2.0, method="shadow_price")
     sm, sc = _shadow_helper(model)
     np.testing.assert_array_equal(pm, sm)
     np.testing.assert_allclose(pc, sc, rtol=1e-6, atol=1e-12)
@@ -477,8 +477,37 @@ def test_dispatcher_falls_back_to_finite_diff_without_usable_duals(caplog, monke
     model.solver = "scipy"
 
     with caplog.at_level(logging.INFO, logger=mod_ccc.__name__):
-        pm, pc = get_conc_control_coeffs(model, fold_change=2.0)
+        pm, pc = get_conc_control_coeffs(model, fold_change=2.0, method="shadow_price")
     fm, fc = _fd_helper(model, fold_change=2.0)
     np.testing.assert_array_equal(pm, fm)
     np.testing.assert_allclose(pc, fc, rtol=1e-6, atol=1e-12)
     assert "returns no usable variable duals" in caplog.text
+
+
+def test_default_method_is_the_finite_difference_and_matches_the_helper():
+    model = _build_toy_ec_model(enzyme_ub=8.0)
+    model.reactions.get_by_id("EX_A").lower_bound = -10.0  # breakpoint setup
+    pm, pc = get_conc_control_coeffs(model, fold_change=2.0)
+    fm, fc = _fd_helper(model, fold_change=2.0)
+    np.testing.assert_array_equal(pm, fm)
+    np.testing.assert_allclose(pc, fc, rtol=1e-12, atol=1e-14)
+
+
+def test_default_method_needs_no_duals(monkeypatch):
+    import importlib
+
+    mod_ccc = importlib.import_module("geckopy.limit_proteins.get_conc_control_coeffs")
+
+    def fail(*args, **kwargs):
+        raise AssertionError("the shadow-price path ran")
+
+    monkeypatch.setattr(mod_ccc, "_shadow_price_coeffs", fail)
+    model = _build_toy_ec_model(enzyme_ub=8.0)
+    enz, coeffs = get_conc_control_coeffs(model)
+    assert enz.any()
+
+
+def test_unknown_method_raises():
+    model = _build_toy_ec_model(enzyme_ub=8.0)
+    with pytest.raises(ValueError, match="method must be"):
+        get_conc_control_coeffs(model, method="dual")
diff --git a/tests/test_review_fixes_r2p.py b/tests/test_review_fixes_r2p.py
index 0faab13..c6ba319 100644
--- a/tests/test_review_fixes_r2p.py
+++ b/tests/test_review_fixes_r2p.py
@@ -59,7 +59,7 @@ def _fixed_reduced_cost(monkeypatch, value):
 def test_a_gain_below_half_the_prediction_gives_zero(solver, monkeypatch):
     model = _with_solver(_build_toy_ec_model(enzyme_ub=5.0), solver)
     _stub_check_solve(monkeypatch, model, 0.3 * 5e-4)
-    enz, coeffs = get_conc_control_coeffs(model)
+    enz, coeffs = get_conc_control_coeffs(model, method="shadow_price")
     assert enz[0]
     assert coeffs[0] == 0.0
 
@@ -68,7 +68,7 @@ def test_a_gain_below_half_the_prediction_gives_zero(solver, monkeypatch):
 def test_a_gain_above_half_the_prediction_keeps_the_reduced_cost(solver, monkeypatch):
     model = _with_solver(_build_toy_ec_model(enzyme_ub=5.0), solver)
     _stub_check_solve(monkeypatch, model, 0.6 * 5e-4)
-    _, coeffs = get_conc_control_coeffs(model)
+    _, coeffs = get_conc_control_coeffs(model, method="shadow_price")
     assert coeffs[0] == pytest.approx(1.0, rel=1e-6)
 
 
@@ -78,7 +78,7 @@ def test_a_gain_below_the_floor_is_noise(solver, monkeypatch):
     model = _with_solver(_build_toy_ec_model(enzyme_ub=5.0), solver)
     _fixed_reduced_cost(monkeypatch, 1e-8)
     _stub_check_solve(monkeypatch, model, 4e-12)
-    _, coeffs = get_conc_control_coeffs(model)
+    _, coeffs = get_conc_control_coeffs(model, method="shadow_price")
     assert coeffs[0] == 0.0
 
 
@@ -87,7 +87,7 @@ def test_a_gain_above_the_floor_is_kept(solver, monkeypatch):
     model = _with_solver(_build_toy_ec_model(enzyme_ub=5.0), solver)
     _fixed_reduced_cost(monkeypatch, 1e-8)
     _stub_check_solve(monkeypatch, model, 6e-12)
-    _, coeffs = get_conc_control_coeffs(model)
+    _, coeffs = get_conc_control_coeffs(model, method="shadow_price")
     assert coeffs[0] == pytest.approx(1e-8)
 
 
@@ -97,7 +97,7 @@ def test_a_reduced_cost_below_1e_10_is_checked_and_kept(solver, monkeypatch):
     model = _with_solver(_build_toy_ec_model(enzyme_ub=500.0, enzyme_coeff_in_R=-2000.0), solver)
     _fixed_reduced_cost(monkeypatch, 5e-11)
     _stub_check_solve(monkeypatch, model, 2.5e-12)
-    enz, coeffs = get_conc_control_coeffs(model)
+    enz, coeffs = get_conc_control_coeffs(model, method="shadow_price")
     assert enz[0]
     assert coeffs[0] == pytest.approx(5e-11)
 
@@ -107,7 +107,7 @@ def test_a_zero_reduced_cost_is_not_checked(solver, monkeypatch):
     model = _with_solver(_build_toy_ec_model(enzyme_ub=5.0), solver)
     _fixed_reduced_cost(monkeypatch, 0.0)
     state = _stub_check_solve(monkeypatch, model, 1.0)
-    _, coeffs = get_conc_control_coeffs(model)
+    _, coeffs = get_conc_control_coeffs(model, method="shadow_price")
     assert coeffs[0] == 0.0
     assert state["calls"] == 1
 
@@ -120,7 +120,7 @@ def test_a_usage_reaction_held_at_its_lower_bound_is_not_checked(solver, monkeyp
     usage.upper_bound = 2.0
     _fixed_reduced_cost(monkeypatch, 1.0)
     state = _stub_check_solve(monkeypatch, model, 1.0)
-    enz, coeffs = get_conc_control_coeffs(model)
+    enz, coeffs = get_conc_control_coeffs(model, method="shadow_price")
     assert enz[0]
     assert coeffs[0] == 0.0
     assert state["calls"] == 1
diff --git a/tests/test_scipy_duals_t3.py b/tests/test_scipy_duals_t3.py
index 463fe16..41c9a72 100644
--- a/tests/test_scipy_duals_t3.py
+++ b/tests/test_scipy_duals_t3.py
@@ -45,7 +45,7 @@ def _flux(model, rxn):
 def test_scipy_reads_the_bound_marginal_and_not_the_finite_difference(monkeypatch):
     model = _scipy_model(enzyme_ub=5.0)
     _forbid_finite_difference(monkeypatch)
-    enz, coeffs = get_conc_control_coeffs(model)
+    enz, coeffs = get_conc_control_coeffs(model, method="shadow_price")
     assert enz[0]
     assert coeffs[0] == pytest.approx(1.0, rel=1e-6)
 
@@ -63,7 +63,7 @@ def test_scipy_gives_zero_for_a_bound_that_does_not_bind_while_the_pool_binds(mo
     rxn = model.reactions.get_by_id("usage_prot_E")
     _forbid_finite_difference(monkeypatch)
     assert _forward_reduced_cost(rxn, _flux(model, rxn)) == 0.0
-    _, coeffs = get_conc_control_coeffs(model)
+    _, coeffs = get_conc_control_coeffs(model, method="shadow_price")
     assert coeffs[0] == 0.0
 
 
@@ -75,7 +75,7 @@ def test_scipy_degenerate_vertex_is_resolved_by_the_verification_step(monkeypatc
     assert _forward_reduced_cost(rxn, _flux(model, rxn)) >= 0.0
     monkeypatch.setattr(mod, "_forward_reduced_cost", lambda rxn, flux=None: 1.0)
     _forbid_finite_difference(monkeypatch)
-    _, coeffs = get_conc_control_coeffs(model)
+    _, coeffs = get_conc_control_coeffs(model, method="shadow_price")
     assert coeffs[0] == 0.0
 
 
@@ -93,8 +93,8 @@ def test_scipy_degenerate_vertex_is_resolved_by_the_verification_step(monkeypatc
 def test_scipy_agrees_with_the_other_solvers(build, other):
     reference = build(other)
     scipy_model = build("scipy")
-    enz_ref, ref = get_conc_control_coeffs(reference)
-    enz_sp, sp = get_conc_control_coeffs(scipy_model)
+    enz_ref, ref = get_conc_control_coeffs(reference, method="shadow_price")
+    enz_sp, sp = get_conc_control_coeffs(scipy_model, method="shadow_price")
     np.testing.assert_array_equal(enz_ref, enz_sp)
     np.testing.assert_allclose(sp, ref, rtol=1e-6, atol=1e-12)
 
@@ -112,7 +112,7 @@ def optimize_without_marginals(self, *args, **kwargs):
     monkeypatch.setattr(problem_class, "optimize", optimize_without_marginals)
     mod._FALLBACK_REPORTED.clear()
     with caplog.at_level(logging.INFO, logger=mod.__name__):
-        pm, pc = get_conc_control_coeffs(model, fold_change=2.0)
+        pm, pc = get_conc_control_coeffs(model, fold_change=2.0, method="shadow_price")
     fm, fc = _finite_difference_coeffs(model, ["E"], 2.0, 0.0)
     np.testing.assert_array_equal(pm, fm)
     np.testing.assert_allclose(pc, fc, rtol=1e-6, atol=1e-12)
@@ -136,6 +136,6 @@ def test_nan_marginal_falls_back(monkeypatch):
     monkeypatch.setattr(mod, "_scipy_upper_marginal", lambda rxn, flux: float("nan"))
     rxn = model.reactions.get_by_id("usage_prot_E")
     assert _forward_reduced_cost(rxn, 5.0) is None
-    enz, coeffs = get_conc_control_coeffs(model)
+    enz, coeffs = get_conc_control_coeffs(model, method="shadow_price")
     _, local = _finite_difference_coeffs(model, ["E"], 2.0, 0.0)
     assert coeffs[0] == pytest.approx(local[0], rel=1e-6)