_mfa_model.py · Scientific Agent Skills

1"""Validated carbon maps and a pinned mfapy steady-state forward model."""

2from __future__ import annotations

3

4from collections import Counter

5import contextlib

6import io

7import json

8import math

9from pathlib import Path

10import re

11

12import numpy as np

13from scipy.linalg import cholesky, solve_triangular

14

15

16class InputError(ValueError):

17 """A scientific input contract was not met."""

18

19

20def require(condition, message):

21 if not condition:

22 raise InputError(message)

23

24

25def fields(value, required, optional=()):

26 require(isinstance(value, dict), "Expected an object")

27 require(set(required) <= value.keys(), f"Missing fields: {set(required) - value.keys()}")

28 require(value.keys() <= set(required) | set(optional),

29 f"Unknown fields: {value.keys() - set(required) - set(optional)}")

30

31

32def number(value, label):

33 require(isinstance(value, (int, float)) and not isinstance(value, bool)

34 and math.isfinite(value), f"{label} must be a finite number")

35 return float(value)

36

37

38def identifier(value):

39 # mfapy generates numerical code internally; only these identifiers reach it.

40 require(isinstance(value, str) and re.fullmatch(r"[A-Za-z][A-Za-z0-9]{0,39}", value),

41 f"Invalid identifier {value!r}; use ASCII letters/digits, starting with a letter")

42 return value

43

44

45def unique(items, label):

46 require(len(set(items)) == len(items), f"Duplicate {label}")

47

48

49def read_json(path):

50 def pairs(items):

51 result = {}

52 for key, value in items:

53 require(key not in result, f"Duplicate JSON key: {key}")

54 result[key] = value

55 return result

56

57 def invalid_constant(value):

58 raise InputError(f"Non-finite JSON constant: {value}")

59

60 return json.loads(Path(path).read_text(encoding="utf-8"), object_pairs_hook=pairs,

61 parse_constant=invalid_constant)

62

63

64class Network:

65 def __init__(self, spec):

66 fields(spec, {"schema_version", "name", "flux_unit", "provenance", "metabolites",

67 "reactions", "fragments"})

68 require(type(spec["schema_version"]) is int and spec["schema_version"] == 1,

69 "Only schema_version 1 is supported")

70 for key in ("name", "flux_unit", "provenance"):

71 require(isinstance(spec[key], str) and bool(spec[key].strip()), f"Missing {key}")

72 for key, maximum in (("metabolites", 100), ("reactions", 200), ("fragments", 100)):

73 require(isinstance(spec[key], list) and 0 < len(spec[key]) <= maximum,

74 f"{key} must be a nonempty list with at most {maximum} entries")

75 self.spec = spec

76 self.metabolites = {}

77 for m in spec["metabolites"]:

78 fields(m, {"id", "carbons", "boundary"}, {"symmetric"})

79 name = identifier(m["id"])

80 require(name not in self.metabolites, f"Duplicate metabolite {name}")

81 require(type(m["carbons"]) is int and 1 <= m["carbons"] <= 12,

82 f"{name}: carbon count must be 1..12")

83 require(m["boundary"] in ("source", "internal", "sink"), f"{name}: invalid boundary")

84 require(type(m.get("symmetric", False)) is bool, f"{name}: symmetric must be boolean")

85 require(not m.get("symmetric") or m["boundary"] == "internal",

86 "Reversal symmetry is supported only for internal pools")

87 self.metabolites[name] = m

88 self.ids = []

89 self.lower, self.upper, self.fixed = [], [], {}

90 balance = {name: [] for name, m in self.metabolites.items() if m["boundary"] == "internal"}

91 self.reactions = {}

92 used = set()

93 for i, r in enumerate(spec["reactions"]):

94 fields(r, {"id", "substrates", "products", "lower", "upper"}, {"fixed"})

95 rid = identifier(r["id"])

96 require(rid not in self.ids, f"Duplicate reaction {rid}")

97 self.ids.append(rid)

98 lo, hi = number(r["lower"], rid), number(r["upper"], rid)

99 require(0 <= lo < hi, f"{rid}: require 0 <= lower < upper; use fixed for fixed fluxes")

100 self.lower.append(lo)

101 self.upper.append(hi)

102 if "fixed" in r:

103 fixed = number(r["fixed"], rid)

104 require(lo <= fixed <= hi, f"{rid}: fixed flux outside bounds")

105 self.fixed[rid] = fixed

106 names, labels = [], []

107 for side in ("substrates", "products"):

108 require(isinstance(r[side], list) and bool(r[side]), f"{rid}: empty {side}")

109 side_names, side_labels = [], []

110 for term in r[side]:

111 require(isinstance(term, list) and len(term) == 2, f"{rid}: expected [metabolite, atoms]")

112 metabolite, atoms = term

113 require(isinstance(metabolite, str) and metabolite in self.metabolites,

114 f"{rid}: unknown metabolite {metabolite!r}")

115 m = self.metabolites[metabolite]

116 require(isinstance(atoms, str) and re.fullmatch(r"[A-Za-z]+", atoms)

117 and len(atoms) == m["carbons"], f"{rid}: atom count mismatch for {metabolite}")

118 require(not (side == "substrates" and m["boundary"] == "sink"),

119 f"{rid}: a sink cannot be consumed")

120 require(not (side == "products" and m["boundary"] == "source"),

121 f"{rid}: a source cannot be produced")

122 side_names.append(metabolite)

123 side_labels.append(atoms)

124 used.add(metabolite)

125 all_atoms = "".join(side_labels)

126 unique(all_atoms, f"carbon atoms on {rid} {side}")

127 names.append(side_names)

128 labels.append(side_labels)

129 require(set("".join(labels[0])) == set("".join(labels[1])),

130 f"{rid}: carbon atoms are lost or created; include CO2 and other carbon products")

131 left, right = Counter(names[0]), Counter(names[1])

132 for name in balance:

133 balance[name].append(right[name] - left[name])

134

135 def stoichiometry(parts):

136 return "+".join((f"{{{count}}}" if count > 1 else "") + name

137 for name, count in Counter(parts).items())

138

139 self.reactions[rid] = dict(

140 order=i, stoichiometry="-->".join(stoichiometry(side) for side in names),

141 reaction="-->".join("+".join(side) for side in names),

142 atommap="-->".join("+".join(side) for side in labels),

143 externalids="", lb=lo, ub=hi)

144 require(used == self.metabolites.keys(), "Unused metabolites must be removed")

145 require(bool(balance), "At least one internal metabolite is required")

146 for name, row in balance.items():

147 require(any(v > 0 for v in row) and any(v < 0 for v in row),

148 f"{name}: internal pool needs production and consumption")

149 self.sources = {name for name, m in self.metabolites.items() if m["boundary"] == "source"}

150 require(bool(self.sources), "At least one carbon source is required")

151 # Catalytic cycles require an initial pool but can receive labeled carbon;

152 # connectivity is checked on the undirected graph below, not reachability.

153 connected = set(self.sources)

154 for _ in self.metabolites:

155 for r in spec["reactions"]:

156 members = {term[0] for term in r["substrates"] + r["products"]}

157 if members & connected:

158 connected.update(members)

159 require(connected == self.metabolites.keys(), "Network contains a pool disconnected from all sources")

160 self.balance = np.array(list(balance.values()), dtype=float)

161 rows, rhs = list(self.balance), [0.] * len(self.balance)

162 for rid, value in self.fixed.items():

163 rows.append(np.eye(len(self.ids))[self.ids.index(rid)])

164 rhs.append(value)

165 self.equalities = np.array(rows)

166 self.rhs = np.array(rhs)

167 self.lower, self.upper = np.array(self.lower), np.array(self.upper)

168 self.fragments = {}

169 for f in spec["fragments"]:

170 fields(f, {"id", "metabolite", "carbons", "provenance"})

171 fid = identifier(f["id"])

172 require(fid not in self.fragments, f"Duplicate fragment {fid}")

173 require(isinstance(f["metabolite"], str) and f["metabolite"] in self.metabolites,

174 f"{fid}: unknown metabolite")

175 m = self.metabolites[f["metabolite"]]

176 require(m["boundary"] == "internal", f"{fid}: fragment must reference an internal pool")

177 require(isinstance(f["carbons"], list) and bool(f["carbons"]), f"{fid}: empty carbon positions")

178 require(all(type(c) is int and 1 <= c <= m["carbons"] for c in f["carbons"]),

179 f"{fid}: carbon positions are 1-based and must exist")

180 unique(f["carbons"], f"positions in {fid}")

181 require(isinstance(f["provenance"], str) and bool(f["provenance"].strip()),

182 f"{fid}: document the measured carbon assignment")

183 self.fragments[fid] = f

184 self._engine = None

185

186 def check_flux(self, flux):

187 require(isinstance(flux, dict) and set(flux) == set(self.ids), "Supply exactly one flux for every reaction")

188 vector = np.array([number(flux[rid], rid) for rid in self.ids])

189 require(self.feasible(vector), "Fluxes violate mass balance, fixed constraints, or bounds")

190 return vector

191

192 def feasible(self, vector, tolerance=1e-8):

193 tolerance *= max(float(np.max(np.abs(vector))), float(np.max(np.abs(self.rhs))), 1e-300)

194 return bool(np.all(np.isfinite(vector)) and np.all(vector >= self.lower - tolerance)

195 and np.all(vector <= self.upper + tolerance)

196 and np.max(np.abs(self.equalities @ vector - self.rhs)) <= tolerance)

197

198 def engine(self):

199 if self._engine is None:

200 import mfapy

201 metabolites = {

202 name: dict(C_number=m["carbons"], order=i,

203 symmetry="symmetry" if m.get("symmetric") else "no",

204 carbonsource="carbonsource" if m["boundary"] == "source" else "no",

205 excreted="excreted" if m["boundary"] == "sink" else "no",

206 externalids="", lb=0., ub=1000.)

207 for i, (name, m) in enumerate(self.metabolites.items())}

208 fragments = {

209 fid: dict(atommap=f["metabolite"] + "_" + ":".join(map(str, sorted(f["carbons"]))),

210 type="intermediate", formula="", use="use", order=i)

211 for i, (fid, f) in enumerate(self.fragments.items())}

212 # The upstream constructor prints diagnostics. Return clean JSON from the CLI.

213 with contextlib.redirect_stdout(io.StringIO()):

214 self._engine = mfapy.metabolicmodel.MetabolicModel(self.reactions, {}, metabolites, fragments)

215 require(hasattr(self._engine, "func"), "mfapy could not construct this carbon network")

216 return self._engine

217

218 def predict(self, vector, experiment):

219 engine = self.engine()

220 # mfapy's generated steady-state solver drops rows with absolute turnover

221 # <= 0.001. MIDs are invariant to uniform flux scaling: keep every positive

222 # reaction at least 1 in the forward solver, without changing fitted units.

223 vector = np.asarray(vector, dtype=float)

224 require(np.all(np.isfinite(vector)) and np.min(vector) >= 0,

225 "Isotope simulation requires finite nonnegative directional fluxes")

226 positive = vector[vector > 0]

227 require(len(positive) > 0, "Isotope steady state is undefined with all fluxes zero")

228 forward = vector / positive.min()

229 cs = engine.generate_carbon_source_template()

230 for name, distribution in experiment["substrates"].items():

231 # Direct arrays avoid upstream set_each_isotopomer's exact sum > 1 check.

232 values = [0.] * (2 ** self.metabolites[name]["carbons"])

233 for bits, fraction in distribution.items():

234 values[int(bits[::-1], 2)] = fraction

235 require(cs.set_all_isotopomers(name, values, correction="no"), f"Invalid source {name}")

236 state = {"reaction": {rid: {"value": value} for rid, value in zip(self.ids, forward)},

237 "metabolite": {name: {"value": 1.} for name in self.metabolites}, "reversible": {}}

238 with np.errstate(all="raise"):

239 mdv = engine.generate_mdv(state, cs)

240 result = {}

241 for fid, f in self.fragments.items():

242 values = np.asarray(mdv.get_fragment_mdv(fid), dtype=float)

243 require(len(values) == len(f["carbons"]) + 1 and np.all(np.isfinite(values))

244 and np.min(values) >= -1e-8 and np.max(values) <= 1 + 1e-8

245 and abs(values.sum() - 1) < 1e-6,

246 f"Invalid simulated MDV for {fid}; check zero-throughput pools and mapping")

247 result[fid] = values

248 return result

249

250

251class Observations:

252 def __init__(self, spec, network):

253 fields(spec, {"schema_version", "steady_state", "shared_fluxes", "mdv_basis",

254 "correction_notes", "experiments"})

255 require(type(spec["schema_version"]) is int and spec["schema_version"] == 1, "Unsupported schema")

256 require(spec["steady_state"] == {"metabolic": True, "isotopic": True}

257 and all(type(v) is bool for v in spec["steady_state"].values()),

258 "This solver requires metabolic AND isotopic steady state; use an INST-MFA engine for time courses")

259 require(spec["shared_fluxes"] is True, "Experiments in a joint fit must share the same flux state")

260 require(spec["mdv_basis"] == "tracer_only", "Correct natural abundance externally; raw MDVs are not accepted")

261 require(isinstance(spec["correction_notes"], str) and bool(spec["correction_notes"].strip()),

262 "Record the correction method and how tracer purity was handled")

263 require(isinstance(spec["experiments"], list) and bool(spec["experiments"]), "No experiments")

264 self.spec, self.network, self.blocks, self.flux_blocks = spec, network, [], []

265 self.warnings = []

266 ids = []

267 for exp in spec["experiments"]:

268 fields(exp, {"id", "substrates", "measurements"}, {"flux_measurements"})

269 ids.append(identifier(exp["id"]))

270 require(isinstance(exp["substrates"], dict) and set(exp["substrates"]) == network.sources,

271 "Every experiment must explicitly describe every source, including unlabeled inputs")

272 for name, distribution in exp["substrates"].items():

273 require(isinstance(distribution, dict) and bool(distribution), f"{name}: missing isotopomers")

274 total = 0.

275 for bits, fraction in distribution.items():

276 require(isinstance(bits, str) and re.fullmatch(r"[01]+", bits)

277 and len(bits) == network.metabolites[name]["carbons"], f"{name}: invalid positional isotopomer")

278 total += number(fraction, name)

279 require(fraction >= 0, f"{name}: negative fraction")

280 require(abs(total - 1) < 1e-10, f"{name}: isotopomer fractions must sum to one; no implicit unlabeled remainder")

281 require(isinstance(exp["measurements"], list), "measurements must be a list")

282 seen = []

283 for obs in exp["measurements"]:

284 fields(obs, {"fragment", "mdv", "omit"}, {"sem", "covariance"})

285 fid = obs["fragment"]

286 require(isinstance(fid, str) and fid in network.fragments, f"Unknown fragment {fid!r}")

287 seen.append(fid)

288 n = len(network.fragments[fid]["carbons"]) + 1

289 require(isinstance(obs["mdv"], list) and len(obs["mdv"]) == n, f"{fid}: expected M+0 through M+{n-1}")

290 y = np.array([number(v, fid) for v in obs["mdv"]])

291 require(np.all(y >= 0) and np.all(y <= 1) and abs(y.sum() - 1) < 1e-6,

292 f"{fid}: MDV must be nonnegative and sum to one; do not silently clip or renormalize")

293 omit = obs["omit"]

294 require(type(omit) is int and 0 <= omit < n, f"{fid}: invalid omitted bin")

295 keep = np.array([i for i in range(n) if i != omit])

296 require(("sem" in obs) != ("covariance" in obs), f"{fid}: supply exactly one of sem or covariance")

297 if "covariance" in obs:

298 cov = np.asarray(obs["covariance"], dtype=float)

299 require(cov.shape == (n, n) and np.all(np.isfinite(cov)), f"{fid}: invalid covariance shape/values")

300 require(np.all(np.diag(cov) > 0), f"{fid}: covariance of retained bins must be positive definite")

301 std = np.sqrt(np.diag(cov))

302 correlation = cov / np.outer(std, std)

303 require(np.max(np.abs(correlation - correlation.T)) <= 1e-10

304 and np.all(np.abs(cov.sum(axis=0)) <= np.sum(np.abs(cov), axis=0) * 1e-10)

305 and np.all(np.abs(cov.sum(axis=1)) <= np.sum(np.abs(cov), axis=1) * 1e-10),

306 f"{fid}: full MDV covariance must be symmetric with zero row/column sums")

307 require(np.linalg.eigvalsh(correlation).min() >= -1e-10,

308 f"{fid}: full MDV covariance must be positive semidefinite")

309 reduced = cov[np.ix_(keep, keep)]

310 else:

311 require(isinstance(obs["sem"], list) and len(obs["sem"]) == n, f"{fid}: invalid sem length")

312 sem = np.array([number(v, fid) for v in obs["sem"]])

313 require(np.all(sem > 0), f"{fid}: sem values must be positive")

314 reduced = np.diag(sem[keep] ** 2)

315 self.warnings.append(f"{exp['id']}/{fid}: diagonal SEM approximation; correlations ignored, intervals and goodness-of-fit approximate")

316 uncertainty = float(np.sqrt(np.diag(reduced)).min())

317 require(abs(y.sum() - 1) <= min(1e-10, uncertainty * 1e-5),

318 f"{fid}: MDV normalization error exceeds numerical/error-scale tolerance; revisit upstream normalization")

319 try:

320 chol = cholesky(reduced, lower=True)

321 except np.linalg.LinAlgError as exc:

322 raise InputError(f"{fid}: covariance of retained bins must be positive definite") from exc

323 self.blocks.append((exp, fid, y, keep, chol))

324 unique(seen, f"fragment measurements in {exp['id']}; combine replicates with an explicit error model")

325 flux_obs = exp.get("flux_measurements", [])

326 require(isinstance(flux_obs, list), "flux_measurements must be a list")

327 flux_ids = []

328 for obs in flux_obs:

329 fields(obs, {"reaction", "value", "sem"})

330 require(obs["reaction"] in network.ids, "Unknown measured reaction")

331 require(number(obs["sem"], "flux sem") > 0, "Flux sem must be positive")

332 number(obs["value"], "measured flux")

333 flux_ids.append(obs["reaction"])

334 self.flux_blocks.append((exp["id"], network.ids.index(obs["reaction"]), obs))

335 unique(flux_ids, f"flux measurements in {exp['id']}")

336 unique(ids, "experiment IDs")

337 self.size = sum(len(block[3]) for block in self.blocks) + len(self.flux_blocks)

338

339 def residuals(self, vector, details=False):

340 predictions = {exp["id"]: self.network.predict(vector, exp) for exp in self.spec["experiments"]}

341 residuals, records = [], []

342 for exp, fid, observed, keep, chol in self.blocks:

343 predicted = predictions[exp["id"]][fid]

344 white = solve_triangular(chol, (predicted - observed)[keep], lower=True)

345 residuals.extend(white)

346 records.append(dict(experiment=exp["id"], fragment=fid, observed=observed.tolist(),

347 predicted=predicted.tolist(), retained_bins=keep.tolist(),

348 whitened_residuals=white.tolist()))

349 for exp_id, index, obs in self.flux_blocks:

350 value = (vector[index] - obs["value"]) / obs["sem"]

351 residuals.append(value)

352 records.append(dict(experiment=exp_id, reaction=obs["reaction"], observed=obs["value"],

353 predicted=float(vector[index]), standardized_residual=float(value)))

354 return (np.array(residuals), records) if details else np.array(residuals)