diff --git a/analysis/sampleCHRR.m b/analysis/sampleCHRR.m
index dfa8f67b..63657764 100644
--- a/analysis/sampleCHRR.m
+++ b/analysis/sampleCHRR.m
@@ -85,10 +85,14 @@
 
 % Particular solution v0 and nullspace basis N: v = v0 + N*x.
 % v0 is ANY particular solution; the sampled flux distribution is invariant to
-% its choice and to the nullspace basis, so MATLAB's '\' (a sparse basic
-% solution) and the Python reference's lstsq (minimum-norm solution) sample the
-% same polytope despite returning different v0.
-v0 = Aeq \ beq;
+% its choice and to the nullspace basis, so the minimum-norm solution used here
+% and the Python reference's lstsq sample the same polytope.
+% lsqminnorm rather than '\': Aeq is square whenever the number of metabolites
+% plus the number of folded-in fixed reactions equals the number of reactions,
+% and it is then rank-deficient, for which '\' returns NaN. Those NaNs travel
+% into the polytope and surface much later as a misleading "empty interior"
+% error from sampleChebyshevCenter.
+v0 = lsqminnorm(Aeq, beq);
 N  = null(Aeq);
 d  = size(N, 2);