Skip to content

rm: run the Arrhenius rules RuleMonkey 3.11.0 misread (#1058) #1221

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
wshlavacek merged 1 commit into from
Oct 23, 2026
Merged

rm: run the Arrhenius rules RuleMonkey 3.11.0 misread (#1058) #1221

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 9 additions & 0 deletions changelog.d/1058.fixed.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
- **`method="rm"` runs Arrhenius rules whose energy patterns carry context
(issue #1058).** RuleMonkey 3.12.0 gives them BNG2's rate, expanding a rule at
load where that is exact and counting `dG` at each firing elsewhere
(richardposner/RuleMonkey#109), so `rm` and `RuleMonkeySession` no longer
refuse an energy pattern with a state, a free or wildcard site, a site named
twice or context beyond one bond, context in the rule's own reactant pattern,
a ring-opening unbinding, or a unidirectional rule named `_reverse_...`.
Shapes RuleMonkey cannot run (a state change, a ring closure, a homodimer,
two bonds, constraints) are still refused, by RuleMonkey, naming the rule.
30 changes: 21 additions & 9 deletions docs/user-guide/network-free.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -200,15 +200,27 @@ reverse half of a reversible rule. Energies are in units of RT: there is no RT
in the formula, and a parameter named `RT` is an ordinary parameter, which a
model in kJ/mol divides by where it uses it (issue #1022).

RuleMonkey expands each energy rule at load, and reads only some energy
patterns as BNG2 does: the bonded pair of sites a rule joins or parts, alone or
with one more bond to a molecule of a third type whose site binds nothing else,
on a rule that names only those two sites. That covers pair and cooperativity
patterns. It drops or misapplies a state, a site that must be free, a wildcard
bond, a site named twice, context further out, and context in the rule itself
(richardposner/RuleMonkey#109), so `method="rm"` refuses such a model by name
until that is fixed upstream (issue #1058). Simulate it with `method="ode"` or
`method="ssa"` on the generated network.
RuleMonkey runs an energy rule that makes one bond between two reactant
patterns, or breaks one within one, between molecules of two types, with no
other operation. It expands the rule at load where that is exact, as for pair
and cooperativity patterns, and everywhere else counts `dG` at each firing as
BNG2's network does per reaction: an energy pattern with a state, a site that
must be free, wildcard bonds, a site named twice or context further out, and
context in the rule's own reactant pattern all run at BNG2's rate
(richardposner/RuleMonkey#109; `method="rm"` refused these until issue #1058).
Counting at each firing costs more per event, and more the larger the complexes
an event touches.

RuleMonkey refuses at load, naming the rule and the reason, an energy rule that
changes a state or adds or deletes a molecule, makes or breaks more than one
bond, closes a ring, binds two molecules of one type, has more than two
reactants, or carries include or exclude constraints, and a rule it cannot
expand whose bond can reach an energy pattern joined only by `.`. Simulate such
a model with `method="ode"` or `method="ssa"` on the generated network. BNG2
writes the reverse half of a reversible rule `X` as a rule named `_reverse_X`,
so two unidirectional rules, `X` and a mirror image named `_reverse_X` with the
same `phi` and `Ea0` parameters, read as one reversible rule, and RuleMonkey
gives the second `1-phi` where BNG2 gives it `phi`.

### Compartments

Expand Down
247 changes: 0 additions & 247 deletions python/bngsim/_network_free_checks.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,15 +24,6 @@
A RuleMonkey session uses :func:`clock_readers` to refuse a segment whose labels
are not the clock the model reads.

Issue #1058: RuleMonkey 3.11.0 gives an Arrhenius rule BNG2's rate,
exp(-(Ea0 + phi*dG)), where every energy pattern the rule's bond is part of is
that bonded pair of sites, alone or with one more bond to a third molecule
whose site binds nothing else, and the rule names only those two sites. It
expands each energy rule at load, and outside that it drops or misapplies the
context (richardposner/RuleMonkey#109), so the rule fires at the wrong rate
with no warning. :func:`refuse_rulemonkey_energy_misreads` refuses such a model
by name until the upstream fix is vendored.

An Arrhenius (energy-based) rate law is refused on NFsim (issue #854's review).
The vendored NFsim expands energy rules at load with three silent defects
(RuleWorld/nfsim#102): it divides every energy by an ``RT`` that BioNetGen does
Expand Down Expand Up @@ -260,244 +251,6 @@ def _children(el: ET.Element, name: str) -> list[ET.Element]:
return [c for c in el if _local(c.tag) == name]


def _graph(el: ET.Element) -> tuple[list[ET.Element], dict[str, tuple[int, str, str]], list]:
"""The molecules of a pattern or species, its sites (component id ->
molecule index, molecule type, component name) and its bonds (pairs of
component ids)."""
molecules = [m for lst in _children(el, "ListOfMolecules") for m in lst]
sites = {}
for i, mol in enumerate(molecules):
for lst in _children(mol, "ListOfComponents"):
for comp in lst:
sites[comp.get("id") or ""] = (i, mol.get("name") or "", comp.get("name") or "")
bonds = [
(b.get("site1") or "", b.get("site2") or "")
for lst in _children(el, "ListOfBonds")
for b in lst
]
return molecules, sites, bonds


def _components(mol: ET.Element) -> list[ET.Element]:
return [c for lst in _children(mol, "ListOfComponents") for c in lst]


def _plain_bond_pattern(pattern: ET.Element) -> bool:
"""Whether a pattern names nothing but explicit bonds between its own
molecules: every component bound once to a partner the pattern lists, with
no state, no site named twice on one molecule, and no compartment."""
molecules, sites, bonds = _graph(pattern)
if not molecules or pattern.get("compartment"):
return False
ends = [site for bond in bonds for site in bond]
for mol in molecules:
comps = _components(mol)
names = [c.get("name") for c in comps]
if mol.get("compartment") or not comps or len(set(names)) != len(names):
return False
for comp in comps:
if comp.get("state") is not None or comp.get("numberOfBonds") != "1":
return False
if ends.count(comp.get("id") or "") != 1:
return False
return all(a in sites and b in sites for a, b in bonds)


Site = tuple[str, str] # (molecule type, component name)


def _bond_partners(root: ET.Element) -> dict[Site, set[Site]]:
"""For every site, the sites a bond in the model can join it to: the bonds
of the seed species and of every rule's patterns. A bond a rule adds, to a
molecule it synthesizes too, is in its product pattern."""
partners: dict[Site, set[Site]] = {}
for el in root.iter():
if _local(el.tag) in ("Species", "ReactantPattern", "ProductPattern"):
_, sites, bonds = _graph(el)
for a, b in bonds:
if a in sites and b in sites:
sa, sb = sites[a][1:], sites[b][1:]
partners.setdefault(sa, set()).add(sb)
partners.setdefault(sb, set()).add(sa)
return partners


def _component_names(root: ET.Element) -> dict[str, list[str]]:
"""Each molecule type's component names, a name repeated as declared."""
names: dict[str, list[str]] = {}
for el in root.iter():
if _local(el.tag) == "MoleculeType":
names[el.get("id") or ""] = [
c.get("id") or "" for lst in _children(el, "ListOfComponentTypes") for c in lst
]
return names


def _arrhenius_centers(root: ET.Element) -> tuple[list[tuple[Site, Site]], list[str], bool]:
"""The bond each Arrhenius rule makes or breaks, as a pair of sites; the
rules RuleMonkey would run on more than those two sites, or with the wrong
products; and whether some Arrhenius rule has a shape RuleMonkey refuses
itself at load (try_expand_arrhenius in third_party/rulemonkey), which then
names it better than this check could."""
centers: list[tuple[Site, Site]] = []
with_context: list[str] = []
rule_names = {r.get("name") for r in root.iter() if _local(r.tag) == "ReactionRule"}
for rule in root.iter():
if _local(rule.tag) != "ReactionRule" or not any(
(rl.get("type") or "").lower() == "arrhenius" for rl in _children(rule, "RateLaw")
):
continue
reactants = [rp for lst in _children(rule, "ListOfReactantPatterns") for rp in lst]
products = [pp for lst in _children(rule, "ListOfProductPatterns") for pp in lst]
molecules: list[ET.Element] = []
sites: dict[str, tuple[int, str, str]] = {}
pattern_of: dict[str, int] = {}
for k, rp in enumerate(reactants):
mols, rp_sites, _ = _graph(rp)
sites.update({c: (len(molecules) + v[0], v[1], v[2]) for c, v in rp_sites.items()})
pattern_of.update(dict.fromkeys(rp_sites, k))
molecules += mols
ops = [op for lst in _children(rule, "ListOfOperations") for op in lst]
# RuleMonkey refuses anything but one AddBond joining two reactant
# patterns, or one DeleteBond inside one, between two molecule types.
if len(ops) != 1 or _local(ops[0].tag) not in ("AddBond", "DeleteBond"):
return [], [], True
a, b = ops[0].get("site1") or "", ops[0].get("site2") or ""
if a not in sites or b not in sites or sites[a][1] == sites[b][1]:
return [], [], True
binding = _local(ops[0].tag) == "AddBond"
if binding and (len(reactants) != 2 or pattern_of[a] == pattern_of[b]):
return [], [], True
if not binding and len(reactants) != 1:
return [], [], True
centers.append((sites[a][1:], sites[b][1:]))
# It builds the expanded rule from the two sites alone: one complex
# from two when binding, two from one when unbinding (a bond inside a
# ring, whose products stay one complex, would part it).
bare = len(products) == (1 if binding else 2) and len(molecules) == 2
bare = bare and not any(rp.get("compartment") for rp in reactants)
for mol in molecules:
comps = _components(mol)
if (
mol.get("compartment")
or len(comps) != 1
or comps[0].get("state") is not None
or (comps[0].get("id") or "") not in (a, b)
):
bare = False
# It takes a rule named _reverse_X for the reverse half of rule X and
# gives it 1-phi, which a one-way rule of that name is not.
name = rule.get("name") or ""
if name.startswith("_reverse_") and name[len("_reverse_") :] not in rule_names:
bare = False
if not bare:
with_context.append(name or rule.get("id") or "?")
return centers, with_context, False


def _expanded_exactly(
pattern: ET.Element,
center: tuple[Site, Site],
partners: dict[Site, set[Site]],
names: dict[str, list[str]],
) -> bool:
"""Whether RuleMonkey's expansion of an energy rule with bond *center*
reads the plain pattern *pattern* as BNG2 does: it is the center's two
sites bonded, with at most one more bond, from one of them to a molecule of
a third type that names only that site, which nothing else binds."""
molecules, sites, bonds = _graph(pattern)
kind = set(center)
on_center = [bd for bd in bonds if {sites[bd[0]][1:], sites[bd[1]][1:]} == kind]
if not on_center:
return True # it holds no bond the rule makes or breaks
if len(on_center) > 1:
return False
core = {sites[on_center[0][0]][0], sites[on_center[0][1]][0]}
others = [bd for bd in bonds if bd != on_center[0]]
leaves = [i for i in range(len(molecules)) if i not in core]
if not leaves:
return not others
if len(leaves) > 1 or len(others) != 1:
return False
a, b = others[0]
leaf_end, core_end = (a, b) if sites[a][0] == leaves[0] else (b, a)
if sites[leaf_end][0] != leaves[0] or sites[core_end][0] not in core:
return False
if sites[leaf_end][1] in {site[0] for site in center}:
return False
if len(_components(molecules[leaves[0]])) != 1:
return False
# RuleMonkey adds one "bound" or "free" copy of the context site to the
# expanded rule, so a type with two sites of its name, or of the center
# site's, matches both variants or counts one embedding for two.
core_type, context = sites[core_end][1], sites[core_end][2]
center_site = next(n for t, n in center if t == core_type)
declared = names.get(core_type, [])
if declared.count(context) != 1 or declared.count(center_site) != 1:
return False
# RuleMonkey tests the context site as bound to anything (energy_expand.cpp).
return partners.get(sites[core_end][1:], set()) <= {sites[leaf_end][1:]}


def rulemonkey_energy_misreads(xml_path: str | Path) -> tuple[list[str], list[str]]:
"""The energy patterns (as BNGL) and the Arrhenius rules (by name) of
*xml_path* that RuleMonkey 3.11.0 expands differently from BNG2 (issue
#1058). Both empty for a model with no Arrhenius rule, for one with an
Arrhenius rule RuleMonkey refuses itself, or when the file cannot be
parsed."""
try:
root = ET.parse(str(xml_path)).getroot()
except (ET.ParseError, OSError):
return [], []
centers, with_context, left_to_rulemonkey = _arrhenius_centers(root)
if left_to_rulemonkey or (not centers and not with_context):
return [], []
partners = _bond_partners(root)
names = _component_names(root)
patterns = []
for ep in root.iter():
if _local(ep.tag) != "EnergyPattern":
continue
graphs = _children(ep, "Pattern") or [ep]
if not all(
_plain_bond_pattern(g)
and all(_expanded_exactly(g, c, partners, names) for c in centers)
for g in graphs
):
patterns.append(ep.get("pattern") or ep.get("id") or "?")
return patterns, with_context


def refuse_rulemonkey_energy_misreads(xml_path: str | Path) -> None:
"""Raise :class:`ModelError` if RuleMonkey would expand an Arrhenius rule of
*xml_path* differently from BNG2 (issue #1058).

RuleMonkey 3.11.0 drops or misapplies an energy pattern's context, and an
Arrhenius rule's own reactant context, when it expands the rule at load
(richardposner/RuleMonkey#109). Delete this check, its calls in
``Simulator`` and ``RuleMonkeySession``, and the helpers only it uses, once
the upstream fix is vendored.
"""
patterns, rules = rulemonkey_energy_misreads(xml_path)
if not patterns and not rules:
return
parts = []
if patterns:
parts.append("energy patterns " + ", ".join(f"'{n}'" for n in patterns))
if rules:
parts.append("rules " + ", ".join(f"'{n}'" for n in rules))
raise ModelError(
"RuleMonkey (method='rm' / 'nf_exact') cannot run this Arrhenius model: when it "
f"expands the energy rules at load it would misread {' and '.join(parts)} "
"(richardposner/RuleMonkey#109), so they would fire at the wrong rate. It reads an "
"energy pattern that is the bonded pair of sites a rule joins or parts, alone or "
"with one more bond to a molecule of a third type whose site binds nothing else, on "
"rules that name only those two sites and join two complexes or part one in two. "
"Simulate it with method='ode' or method='ssa' on the generated network "
"(issue #1058)."
)


def compartments(xml_path: str | Path) -> list[str]:
"""The ids of the compartments *xml_path* declares. Empty when the file
cannot be parsed."""
Expand Down
2 changes: 0 additions & 2 deletions python/bngsim/_rulemonkey_session.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -101,12 +101,10 @@ def __init__(

from bngsim._network_free_checks import (
refuse_rulemonkey_compartments,
refuse_rulemonkey_energy_misreads,
rulemonkey_clock,
)

refuse_rulemonkey_compartments(self._xml_path)
refuse_rulemonkey_energy_misreads(self._xml_path) # until issue #1058
# The functions that read the clock, for which a segment's labels must
# be the session clock (see simulate). Where the model reads it through
# its functions alone, the core reads it as time() + __bngsim_t0, set
Expand Down
2 changes: 0 additions & 2 deletions python/bngsim/_simulator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1045,13 +1045,11 @@ def __init__(
from bngsim._bngsim_core import RuleMonkeySimulator
from bngsim._network_free_checks import (
refuse_rulemonkey_compartments,
refuse_rulemonkey_energy_misreads,
rulemonkey_clock,
)
from bngsim._rulemonkey_session import warn_rulemonkey_diagnostics

refuse_rulemonkey_compartments(xml_path)
refuse_rulemonkey_energy_misreads(xml_path) # until issue #1058
# Elapsed time, which keeps event times small at a large t_start. A
# model that reads the clock reads model time (issue #1022): through
# time(), from a copy of its XML whose clock is time() + __bngsim_t0,
Expand Down
Loading
Loading