Examples#
Complete, runnable scripts from
py-bridge/python/examples. Each is a standalone
program — run it with the project’s Python interpreter:
python examples/multiphase_serial_passage.py
Each script’s module docstring explains the chemistry and the modeling choices; read it alongside the code.
Serial passage of a self-replicating species#
Multi-phase stochastic simulation with MultiPhaseStochsim:
a species replicates in a flask phase and is periodically diluted into a
waste phase, producing the characteristic passage “sawtooth”. Shows discrete,
operator-scheduled transfer moves between simulate segments and the exact
mass-balance invariants they satisfy.
1"""Worked example: serial passage of a self-replicating species (MultiPhaseStochsim).
2
3Serial passage is the classic microbiology protocol -- grow a culture, transfer a
4small aliquot into fresh medium, repeat -- used to propagate, adapt, or attenuate a
5population over many generations. Here a species ``A`` replicates (``A -> 2A``) in a
6``flask`` phase; after each growth interval we *dilute* by discarding most of the
7population into a ``waste`` phase, keeping only a 1:``dilution`` aliquot. That
8discard-and-keep step is exactly ``MultiPhaseStochsim.transfer`` with a binomial
9``fraction=`` -- a discrete, operator-scheduled move applied between ``simulate``
10segments, not an intrinsic reaction.
11
12Choosing the growth time so that one interval of growth undoes one dilution
13(``e^{r*T} = dilution``) yields the characteristic passage "sawtooth": the
14population climbs during growth, is cut back at each transfer, and climbs again from
15the survivors.
16
17The only assertions here are the exact invariants a transfer must satisfy -- it
18neither creates nor destroys molecules (mass balance), and what stays plus what is
19discarded equals what was grown. The population sizes themselves are stochastic and
20simply printed.
21
22Run: python examples/multiphase_serial_passage.py
23"""
24import math
25import os
26import sys
27
28# Run standalone from any working directory: put the molfoundry package (one
29# directory up, in py-bridge/python) on the import path.
30sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
31
32from molfoundry import MultiPhaseStochsim, graphGMLString, ruleGMLString # noqa: E402
33
34
35def _species(label: str):
36 return graphGMLString(
37 'graph [ node [ id 0 label "%s" ] ]' % label, label, add=False
38 )
39
40
41def _duplicate(label: str):
42 # A -> A + A: a second copy appears on the right-hand side (autocatalytic
43 # growth). Its mass-action propensity is rate * n, so the population grows
44 # exponentially at `rate` per capita.
45 return ruleGMLString(
46 'rule [ ruleID "dup_%s" '
47 'left [ node [ id 0 label "%s" ] ] '
48 'right [ node [ id 0 label "%s" ] node [ id 1 label "%s" ] ] ]'
49 % (label, label, label, label),
50 "%s->2%s" % (label, label),
51 add=False,
52 )
53
54
55def main() -> None:
56 a = _species("A")
57 growth_rate = 1.0
58 dilution = 10 # 1:10 passage -> keep fraction 1/10
59 keep = 1.0 / dilution
60 # Grow just enough to undo the dilution each cycle: e^{r*T} = dilution.
61 growth_time = math.log(dilution) / growth_rate
62 passages = 6
63 seed = 20260808
64
65 mp = MultiPhaseStochsim(seed=seed)
66 flask = mp.addPhase(
67 "flask",
68 graphDatabase=[a],
69 expandStrategy=_duplicate("A"),
70 initialState={a: 50},
71 reactionRate=growth_rate,
72 )
73 # A sink for the majority poured off at each passage (models discarding it).
74 waste = mp.addPhase(
75 "waste", graphDatabase=[a], expandStrategy=[], initialState={a: 0}
76 )
77
78 print(
79 "serial passage: 1:%d dilution, growth_time=%.3f, seed=%d\n"
80 % (dilution, growth_time, seed)
81 )
82 print("%7s %7s %7s %9s %6s" % ("passage", "grown", "aliquot", "discarded", "ratio"))
83 for p in range(1, passages + 1):
84 mp.simulate(time=growth_time) # grow A -> 2A for one interval
85 grown = mp.state(flask, a)
86 before = mp.total(a) # flask + waste, for the mass-balance check
87 discarded = mp.transfer(flask, waste, species=a, fraction=1.0 - keep)
88 aliquot = mp.state(flask, a)
89 # A transfer moves molecules; it never creates or destroys them.
90 assert mp.total(a) == before
91 assert aliquot + discarded == grown
92 ratio = aliquot / grown if grown else 0.0
93 print("%7d %7d %7d %9d %6.3f" % (p, grown, aliquot, discarded, ratio))
94
95 print("\nkept-fraction target = %.3f (1:%d)" % (keep, dilution))
96 print(
97 "total A ever (flask + waste) = %d "
98 "-- grew by replication, conserved across every transfer" % mp.total(a)
99 )
100
101
102if __name__ == "__main__":
103 main()
Enzyme mechanism → electron-flow rules#
Translating a curated M-CSA enzyme mechanism into molfoundry electron-graph
rules (Rule over GraphType.Electron). Each curly arrow of
the mechanism becomes one graph rewrite that relocates electron ownership; the
example parses the ChemAxon MEFlow arrows and emits composable per-step rules.
1"""Translate an M-CSA enzyme mechanism into molfoundry electron-graph rules.
2
3The Mechanism and Catalytic Site Atlas (M-CSA, https://www.ebi.ac.uk/thornton-srv/m-csa/)
4publishes curated enzyme reaction mechanisms as a sequence of steps, each drawn
5as a 2D scheme with electron-flow ("curly") arrows. Each step is downloadable as
6a ChemAxon Marvin document (``.mrv``) whose ``<MEFlow>`` elements encode the
7arrows: the tail is a lone pair or a bond, the head an atom or a bond.
8
9That formalism maps exactly onto molfoundry's half-edge electron graph
10(``GraphType.Electron``), where every electron is a vertex labeled ``<e>``, an
11``atom--<e>`` edge is ownership (the sigma map) and an ``<e>--<e>`` edge a shared
12pair. A curly arrow is one graph rewrite: relocate the ownership of one electron
13of a pair from the losing atom to the gaining atom -- delete the old ``atom--<e>``
14edge in ``left``, create the new one in ``right``, keep everything else in
15``context``. Every M-CSA arrow is one of three cases:
16
17 A lone pair X -> bond X-Y (X donates: e_move X -> Y)
18 B bond X-Y -> lone pair Z (collapse onto Z: e_move other -> Z)
19 C bond X-Y (pivot X) -> bond X-Z (pi shift: e_move Y -> Z)
20
21An arrow only says which electrons move, so a bare per-step rule is an
22under-constrained motif. Each rule therefore carries a selectable amount of
23unchanged scaffold as context (see :class:`RuleContext`): the default keeps the
24unchanged bonds among the reacting atoms (so, e.g., a pi-shift promotes an
25existing single bond rather than bonding two unbonded atoms), and higher levels
26add a neighbor shell or the whole active-site skeleton.
27
28This example downloads one entry (glutamate racemase, M-CSA #1), converts every
29mechanism step into one electron-graph rule, composes the per-step rules into a
30single "enzyme" rule with molfoundry's rule-composition engine, and prints each
31rule and the composite with ``GraphPrinter`` in circles mode. The rendered SVG
32summary is written to ``./summary/summary.html``.
33
34Run (from py-bridge/python, with the dev venv active):
35
36 python examples/mcsa_electron_rules.py # entry 1, PARTICIPANTS context
37 python examples/mcsa_electron_rules.py 2 # any M-CSA entry id
38 python examples/mcsa_electron_rules.py 1 2 # context level: 0=MOTIF .. 3=SCENE
39
40The pure parsing/translation helpers (``parse_mrv``, ``classify_arrow``,
41``build_global_ids``, ``emit_electron_rule_gml`` ...) are framework-free and are
42unit-tested in tests/test_mcsa_electron_rules.py.
43"""
44import enum
45import json
46import os
47import re
48import sys
49import urllib.request
50import xml.etree.ElementTree as ET
51from collections import Counter
52
53# Run standalone from any working directory: put the molfoundry package (one
54# directory up, in py-bridge/python) on the import path.
55sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
56
57from molfoundry import ( # noqa: E402
58 GraphType,
59 RCEvaluator,
60 enable_post,
61 postSection,
62 rcSuper,
63 ruleGMLString,
64 smiles,
65)
66from molfoundry.graphprinter import GraphPrinter # noqa: E402
67
68MCSA_API = "https://www.ebi.ac.uk/thornton-srv/m-csa/api/entries/%s/?format=json"
69CACHE_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), ".mcsa_cache")
70MRV_NS = "{http://www.chemaxon.com}"
71
72
73class RuleContext(enum.IntEnum):
74 """How much unchanged scaffold a per-step rule carries around its arrows.
75
76 An arrow only defines *which electrons move*; on its own that is an
77 under-constrained rewrite (e.g. a pi-shift with no sigma bond would form a
78 bond between unbonded atoms). Adding the unchanged bonds as context anchors
79 the rewrite and makes it more chemically specific.
80 """
81
82 MOTIF = 0 # only the pushed electrons -- a bare mechanistic motif
83 PARTICIPANTS = 1 # + unchanged bonds among the arrow-touched atoms (default)
84 ENVIRONMENT = 2 # + one shell of neighboring atoms and their bonds
85 SCENE = 3 # + the whole drawn active-site skeleton
86
87
88# ---------------------------------------------------------------------------
89# Download + cache (the only impure part of the loader)
90# ---------------------------------------------------------------------------
91def _get(url, binary=False):
92 req = urllib.request.Request(url, headers={"User-Agent": "molfoundry-example"})
93 with urllib.request.urlopen(req, timeout=30) as resp:
94 data = resp.read()
95 return data if binary else data.decode("utf-8")
96
97
98def load_entry(mcsa_id, cache_dir=CACHE_DIR):
99 """Fetch one M-CSA entry and its per-step Marvin documents.
100
101 Everything is cached under ``cache_dir`` so re-runs (and offline runs) do no
102 network I/O. Returns ``(enzyme_name, [(step_id, description, is_product,
103 mrv_text), ...])`` for the first detailed mechanism.
104 """
105 entry_dir = os.path.join(cache_dir, "mcsa_%s" % mcsa_id)
106 os.makedirs(entry_dir, exist_ok=True)
107
108 entry_path = os.path.join(entry_dir, "entry.json")
109 if os.path.exists(entry_path):
110 with open(entry_path, "r", encoding="utf-8") as fh:
111 entry = json.load(fh)
112 else:
113 entry = json.loads(_get(MCSA_API % mcsa_id))
114 with open(entry_path, "w", encoding="utf-8") as fh:
115 json.dump(entry, fh)
116
117 mechanisms = entry["reaction"]["mechanisms"]
118 mech = next((m for m in mechanisms if m.get("is_detailed")), mechanisms[0])
119
120 steps = []
121 for step in mech["steps"]:
122 url = step["marvin_xml"]
123 if not url.startswith("http"):
124 url = "https://" + url
125 mrv_path = os.path.join(entry_dir, "step_%s.mrv" % step["step_id"])
126 if os.path.exists(mrv_path):
127 with open(mrv_path, "r", encoding="utf-8") as fh:
128 mrv = fh.read()
129 else:
130 mrv = _get(url)
131 with open(mrv_path, "w", encoding="utf-8") as fh:
132 fh.write(mrv)
133 steps.append((step["step_id"], step["description"],
134 bool(step.get("is_product")), mrv))
135 return entry["enzyme_name"], steps
136
137
138# ---------------------------------------------------------------------------
139# Pure MRV parsing
140# ---------------------------------------------------------------------------
141def parse_mrv(text):
142 """Parse one Marvin step document.
143
144 Returns a dict with ``atoms`` (id -> {el, x, y, fc}), ``bonds``
145 (frozenset(id, id) -> order) and ``arrows`` (list of ``[src_ids, dst_ids]``,
146 where each id list has one element for a lone pair / atom and two for a bond).
147 """
148 root = ET.fromstring(text)
149 atoms = {}
150 for atom in root.iter(MRV_NS + "atom"):
151 atoms[atom.get("id")] = {
152 "el": atom.get("elementType"),
153 "x": round(float(atom.get("x2")), 1),
154 "y": round(float(atom.get("y2")), 1),
155 "fc": int(atom.get("formalCharge") or 0),
156 }
157 bonds = {}
158 for bond in root.iter(MRV_NS + "bond"):
159 i, j = bond.get("atomRefs2").split()
160 bonds[frozenset((i, j))] = int(bond.get("order") or 1)
161 arrows = []
162 for flow in root.iter(MRV_NS + "MEFlow"):
163 points = []
164 for child in flow:
165 refs = child.get("atomRefs") or child.get("atomRef") or ""
166 points.append([r.split(".")[-1] for r in refs.split()])
167 arrows.append(points)
168 return {"atoms": atoms, "bonds": bonds, "arrows": arrows}
169
170
171def classify_arrow(arrow):
172 """Resolve one MEFlow arrow to ``(kind, keep, from_atom, to_atom)``.
173
174 Every arrow moves one electron of a pair from ``from_atom`` to ``to_atom``;
175 ``keep`` owns the pair's stationary electron. See the three cases A/B/C in the
176 module docstring.
177 """
178 src, dst = arrow[0], arrow[1]
179 if len(src) == 1 and len(dst) == 2: # A: lone pair -> bond
180 x = src[0]
181 y = dst[0] if dst[1] == x else dst[1]
182 return "A", x, x, y
183 if len(src) == 2 and len(dst) == 1: # B: bond -> lone pair
184 z = dst[0]
185 other = src[0] if src[1] == z else src[1]
186 return "B", z, other, z
187 if len(src) == 2 and len(dst) == 2: # C: bond -> bond (pi shift)
188 pivot = (set(src) & set(dst)).pop()
189 y = src[0] if src[1] == pivot else src[1]
190 z = dst[0] if dst[1] == pivot else dst[1]
191 return "C", pivot, y, z
192 raise ValueError("unclassifiable MEFlow arrow: %r" % arrow)
193
194
195# ---------------------------------------------------------------------------
196# Global atom identity across steps
197# ---------------------------------------------------------------------------
198def build_global_ids(step_atoms):
199 """Assign every atom a stable global id shared across steps.
200
201 M-CSA draws each step from the same layout, so a heavy atom keeps its 2D
202 coordinate between steps; only the atom(s) that move (a transferred proton, an
203 occasionally redrawn charged center) change position, and steps that renumber
204 are re-registered by coordinate. Atoms are matched step-to-step by
205 (element, coordinate) against a running registry that follows the movers;
206 an atom that moved this step is reconciled to the leftover global atom of its
207 element. Returns a list (one per step) of ``{local_id: global_id}`` maps.
208 """
209 registry = {} # (el, x, y) -> global_id, refreshed to current positions
210 gid_el = {} # global_id -> element
211 pos = {} # global_id -> (x, y), current position
212 maps = []
213 next_gid = 0
214 for atoms in step_atoms:
215 local_map, used = {}, set()
216 unmatched = []
217 for lid, a in atoms.items():
218 gid = registry.get((a["el"], a["x"], a["y"]))
219 if gid is not None and gid not in used:
220 local_map[lid] = gid
221 used.add(gid)
222 else:
223 unmatched.append(lid)
224 # An unmatched atom either moved (reuse a leftover global of same element)
225 # or is genuinely new (fresh global id).
226 leftover = [g for g in pos if g not in used]
227 for lid in unmatched:
228 el = atoms[lid]["el"]
229 same = [g for g in leftover if gid_el[g] == el]
230 if same:
231 gid = same[0]
232 leftover.remove(gid)
233 else:
234 gid = next_gid
235 next_gid += 1
236 gid_el[gid] = el
237 local_map[lid] = gid
238 used.add(gid)
239 # Refresh the registry to this step's positions so the next step matches
240 # movers at their new location.
241 registry = {}
242 for lid, gid in local_map.items():
243 a = atoms[lid]
244 pos[gid] = (a["x"], a["y"])
245 registry[(a["el"], a["x"], a["y"])] = gid
246 maps.append(local_map)
247 return maps, gid_el
248
249
250# ---------------------------------------------------------------------------
251# Chemical self-check: applying a step's arrows must reproduce the next step
252# ---------------------------------------------------------------------------
253def apply_arrows_chemical(parsed):
254 """Apply a step's arrows to its drawn structure, returning (bonds, charges).
255
256 Pure electron bookkeeping on the chemical graph -- used only to validate the
257 translation against M-CSA's own next-step drawing.
258 """
259 bonds = dict(parsed["bonds"])
260 charge = {aid: parsed["atoms"][aid]["fc"] for aid in parsed["atoms"]}
261
262 def bump(i, j, delta):
263 k = frozenset((i, j))
264 bonds[k] = bonds.get(k, 0) + delta
265 if bonds[k] == 0:
266 del bonds[k]
267
268 for arrow in parsed["arrows"]:
269 kind, keep, frm, to = classify_arrow(arrow)
270 if kind == "A": # lone pair on `keep` forms keep-to bond
271 bump(keep, to, +1)
272 charge[keep] += 1
273 charge[to] -= 1
274 elif kind == "B": # frm-keep bond collapses onto `keep`
275 bump(frm, keep, -1)
276 charge[frm] += 1
277 charge[keep] -= 1
278 else: # C: keep-frm bond shifts to keep-to
279 bump(keep, frm, -1)
280 bump(keep, to, +1)
281 charge[frm] += 1
282 charge[to] -= 1
283 return bonds, charge
284
285
286def chemical_signature(atoms, bonds, charge):
287 """Order-independent multiset signature of a chemical state."""
288 degree = Counter()
289 for edge, order in bonds.items():
290 for aid in edge:
291 degree[aid] += order
292 atom_sig = Counter((atoms[a]["el"], charge[a], degree[a]) for a in atoms)
293 bond_sig = Counter()
294 for edge, order in bonds.items():
295 i, j = tuple(edge)
296 bond_sig[(tuple(sorted((atoms[i]["el"], atoms[j]["el"]))), order)] += 1
297 return atom_sig, bond_sig
298
299
300def self_check(parsed_steps):
301 """Verify each non-product step's arrows reproduce the next step's chemistry.
302
303 A step whose successor drops spectator fragments (e.g. the released product)
304 is only checked on the atoms the two steps share.
305 """
306 results = []
307 for i in range(len(parsed_steps) - 1):
308 cur, nxt = parsed_steps[i], parsed_steps[i + 1]
309 if not cur["arrows"]:
310 continue
311 bonds, charge = apply_arrows_chemical(cur)
312 shared = set(cur["atoms"]) & set(nxt["atoms"])
313 # Compare only over atoms present in both drawings, and only when the two
314 # steps actually share their id space (no renumbering between them).
315 renumbered = any(cur["atoms"][a]["el"] != nxt["atoms"][a]["el"] for a in shared)
316 if renumbered or len(shared) < len(cur["atoms"]):
317 results.append((i + 1, None)) # not directly comparable
318 continue
319 got = chemical_signature(cur["atoms"], bonds, charge)
320 want_bonds = nxt["bonds"]
321 want_charge = {a: nxt["atoms"][a]["fc"] for a in nxt["atoms"]}
322 want = chemical_signature(nxt["atoms"], want_bonds, want_charge)
323 results.append((i + 1, got == want))
324 return results
325
326
327# ---------------------------------------------------------------------------
328# Electron-graph rule emission
329# ---------------------------------------------------------------------------
330def step_moves(parsed, local_to_global):
331 """The list of ``(kind, keep, from, to)`` moves for a step, in global ids."""
332 moves = []
333 for arrow in parsed["arrows"]:
334 kind, keep, frm, to = classify_arrow(arrow)
335 moves.append((kind, local_to_global[keep], local_to_global[frm],
336 local_to_global[to]))
337 return moves
338
339
340def emit_electron_rule_gml(rule_id, moves, node_id, label_of,
341 context=RuleContext.PARTICIPANTS, bonds=None, atoms=None):
342 """Build the GML of one electron-graph rule from a step's electron moves.
343
344 ``node_id`` maps a global atom id to the integer node id used in the rule and
345 ``label_of`` to its vertex label (element symbol, or a unique tag for
346 composition). ``context`` (see :class:`RuleContext`) selects how much
347 unchanged scaffold to carry as context pairs: ``bonds`` (global
348 frozenset->order) is required above ``MOTIF``, ``atoms`` (global->element)
349 for ``SCENE``. Electron vertices get ids in a private high range.
350 """
351 if context != RuleContext.MOTIF and bonds is None:
352 raise ValueError("bonds are required for context above MOTIF")
353 if context == RuleContext.SCENE and atoms is None:
354 raise ValueError("atoms are required for SCENE context")
355
356 participants = {a for _, k, f, t in moves for a in (k, f, t)}
357 if context <= RuleContext.PARTICIPANTS:
358 included = set(participants)
359 elif context == RuleContext.ENVIRONMENT: # + one neighbor shell
360 included = set(participants) | {v for e in bonds for v in e if e & participants}
361 else: # SCENE: everything drawn
362 included = set(atoms)
363
364 # A bond order not consumed by a B/C arrow (whose source bond is (from, keep))
365 # survives the step and is carried as that many unchanged context pairs.
366 consumed = Counter(frozenset((f, k)) for kind, k, f, t in moves if kind in "BC")
367
368 ctx_nodes, ctx_edges, left_edges, right_edges = [], [], [], []
369 for g in sorted(included, key=node_id):
370 ctx_nodes.append('node [ id %d label "%s" ]' % (node_id(g), label_of(g)))
371
372 electron = 1000
373 for _, keep, frm, to in moves: # the moving pairs (L->R)
374 e_keep, e_move = electron, electron + 1
375 electron += 2
376 ctx_nodes += ['node [ id %d label "<e>" ]' % e_keep,
377 'node [ id %d label "<e>" ]' % e_move]
378 ctx_edges.append("edge [ source %d target %d ]" % (node_id(keep), e_keep))
379 ctx_edges.append("edge [ source %d target %d ]" % (e_keep, e_move))
380 left_edges.append("edge [ source %d target %d ]" % (node_id(frm), e_move))
381 right_edges.append("edge [ source %d target %d ]" % (node_id(to), e_move))
382
383 if context != RuleContext.MOTIF: # unchanged context bonds
384 for edge, order in bonds.items():
385 if not edge <= included:
386 continue
387 u, v = tuple(edge)
388 for _ in range(order - consumed.get(edge, 0)):
389 e_keep, e_move = electron, electron + 1
390 electron += 2
391 ctx_nodes += ['node [ id %d label "<e>" ]' % e_keep,
392 'node [ id %d label "<e>" ]' % e_move]
393 ctx_edges += [
394 "edge [ source %d target %d ]" % (node_id(u), e_keep),
395 "edge [ source %d target %d ]" % (e_keep, e_move),
396 "edge [ source %d target %d ]" % (node_id(v), e_move),
397 ]
398
399 return (
400 'rule [\n ruleID "%s" type "electron"\n'
401 " left [ %s ]\n"
402 " context [ %s ]\n"
403 " right [ %s ]\n]"
404 % (rule_id, " ".join(left_edges),
405 " ".join(ctx_nodes + ctx_edges), " ".join(right_edges))
406 )
407
408
409def globalize_step(parsed, local_to_global):
410 """A step's bonds and atoms re-keyed to global ids: (bonds, atoms)."""
411 bonds = {frozenset((local_to_global[u], local_to_global[v])): order
412 for edge, order in parsed["bonds"].items() for u, v in [tuple(edge)]}
413 atoms = {local_to_global[lid]: a["el"] for lid, a in parsed["atoms"].items()}
414 return bonds, atoms
415
416
417def relabel_to_elements(rule):
418 """Strip the unique ``El_<n>`` composition tags back to plain element labels."""
419 gml = re.sub(r'label "([A-Z][a-z]?)_\d+"', r'label "\1"', rule.getGMLString())
420 return ruleGMLString(gml, name=rule.name, add=False)
421
422
423def compose_enzyme_rule(moves_per_step, node_id, gid_el):
424 """Compose the per-step rules into one rule with molfoundry's RCEvaluator.
425
426 Composition runs on the ``MOTIF`` core (the pushed electrons only): the
427 mechanism's pushes chain, while the surrounding scaffold is unchanged context
428 that only over-constrains the partial-overlap search (carrying it makes
429 ``rcSuper`` fail to chain the steps). The net transformation is the same
430 either way -- context pairs never enter ``left``/``right``.
431
432 Electron vertices share the label ``<e>``, so a naive composition of
433 element-labeled rules is ambiguous (it may glue two distinct sulfurs). We
434 compose over uniquely-labeled copies -- each atom tagged with its global
435 identity -- which forces the one chemically-correct gluing, then relabel the
436 result back to element symbols.
437 """
438 unique_label = lambda g: "%s_%d" % (gid_el[g], node_id(g))
439 rules = [
440 ruleGMLString(emit_electron_rule_gml("step%d" % (i + 1), moves, node_id,
441 unique_label, context=RuleContext.MOTIF),
442 add=False)
443 for i, moves in enumerate(moves_per_step)
444 ]
445 working = [rules[0]]
446 for nxt in rules[1:]:
447 candidates, seen, unique = [], set(), []
448 for w in working:
449 candidates.extend(RCEvaluator().eval(w * rcSuper * nxt))
450 for c in candidates:
451 key = c._ref.getCanonicalKey()
452 if key not in seen:
453 seen.add(key)
454 unique.append(c)
455 if not unique:
456 raise RuntimeError("rule composition produced no result")
457 working = unique
458 # A correct sequential mechanism collapses to a single composite.
459 composite = min(working, key=lambda r: len(re.findall(r"node \[", r.getGMLString())))
460 return relabel_to_elements(composite)
461
462
463# ---------------------------------------------------------------------------
464# Main
465# ---------------------------------------------------------------------------
466def main(mcsa_id, context=RuleContext.PARTICIPANTS):
467 enzyme, raw_steps = load_entry(mcsa_id)
468 print("M-CSA #%s: %s -- %d mechanism steps [rule context: %s]" % (mcsa_id, enzyme, len(raw_steps), context.name))
469
470 parsed_steps = [parse_mrv(mrv) for (_, _, _, mrv) in raw_steps]
471
472 # Validate the arrow translation against M-CSA's own next-step drawings.
473 print("\nSelf-check (arrows reproduce the next step's chemistry):")
474 for step_no, ok in self_check(parsed_steps):
475 label = "n/a (spectators differ)" if ok is None else ("PASS" if ok else "FAIL")
476 print(" step %d -> %d: %s" % (step_no, step_no + 1, label))
477 if ok is False:
478 raise SystemExit("arrow translation is inconsistent with M-CSA")
479
480 # Global atom identity, then one electron-rule per reactive step (product steps
481 # carry no arrows).
482 id_maps, gid_el = build_global_ids([p["atoms"] for p in parsed_steps])
483 reactive_idx = [i for i, p in enumerate(parsed_steps) if p["arrows"]]
484 moves_per_step = [step_moves(parsed_steps[i], id_maps[i]) for i in reactive_idx]
485 globals_per_step = [globalize_step(parsed_steps[i], id_maps[i]) for i in reactive_idx]
486
487 # node_id must cover every atom any context level can pull in, not just the reactive core.
488 node_id = {g: i for i, g in enumerate(sorted({g for m in id_maps for g in m.values()}))}.__getitem__
489 core = sorted({a for moves in moves_per_step for _, k, f, t in moves for a in (k, f, t)}, key=node_id)
490 print("\nReactive core: %d atoms (%s)" % (len(core), ", ".join(gid_el[g] for g in core)))
491
492 enable_post()
493 printer = GraphPrinter()
494 printer.drawMode = "circles"
495
496 # The substrate as an electron graph, in circles mode, for orientation.
497 postSection("Substrate (L-glutamate) as an electron graph")
498 smiles("OC(=O)C(N)CCC(=O)O", name="L-glutamate", graphType=GraphType.Electron).print(printer)
499
500 # One electron-graph rule per mechanism step, at the requested context level.
501 for i, moves, (gbonds, gatoms) in zip(reactive_idx, moves_per_step, globals_per_step):
502 sid, desc = raw_steps[i][0], raw_steps[i][1]
503 gml = emit_electron_rule_gml("step %d" % sid, moves, node_id, lambda g: gid_el[g], context=context,
504 bonds=gbonds, atoms=gatoms)
505 rule = ruleGMLString(gml, name="Step %d" % sid, add=False)
506 postSection("Step %d: %s" % (sid, desc))
507 rule.print(printer)
508
509 # Compose every step into a single "enzyme" rule (on the electron-push core).
510 print("\nComposing %d step rules into one enzyme rule..." % len(moves_per_step))
511 enzyme_rule = compose_enzyme_rule(moves_per_step, node_id, gid_el)
512 postSection("Composed enzyme rule (%s)" % enzyme)
513 enzyme_rule.print(printer)
514 print("Composed rule: type=%s, %d electron relocations net" % (enzyme_rule.type, enzyme_rule.numLeftComponents))
515
516 print("\nWrote SVG summary to ./summary/summary.html")
517
518
519if __name__ == "__main__":
520 entry_id = sys.argv[1] if len(sys.argv) > 1 else "1"
521 ctx = RuleContext(int(sys.argv[2])) if len(sys.argv) > 2 else RuleContext.SCENE
522 main(entry_id, ctx)