From 8458f1523acad3b2be11c9adc314f3babe6251d2 Mon Sep 17 00:00:00 2001 From: Joel Eliason Date: Mon, 15 Jun 2026 15:38:02 -0600 Subject: [PATCH] Support on constant=false parameters (promote to state) SimBiology emits an abstract non-physical state (e.g. a dimensionless niche fraction) as a parameter with ConstantValue=false plus a rate rule, rather than as a species. The generator only parsed , so a targeting a parameter was silently dropped: the parameter stayed frozen at its initial value with no derivative emitted, and any downstream coupling to it evaluated against a constant. _parse_rate_rules promotes each such parameter to an amount-tracked pseudo-species (state slot + SP_ enum + initial value) and injects one synthetic production reaction whose rate law IS the rate-rule expression (d(var)/dt = expr). The existing stoichiometry -> ydot -> Jacobian -> RHS machinery then evolves it like any other state, with no special-casing downstream. Rate rules whose target is a species (already a state) are left to the normal reaction path. Verified on a PDAC model whose TLA niche fraction is the model's only rateRule: the state now integrates and its H_CXCL13-driven formation + decay are wired into both the RHS and the analytical Jacobian. --- src/qsp_codegen/codegen.py | 50 +++++++++++++++++++++ tests/test_rate_rules.py | 92 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 142 insertions(+) create mode 100644 tests/test_rate_rules.py diff --git a/src/qsp_codegen/codegen.py b/src/qsp_codegen/codegen.py index 9dcef2f..4d1c96c 100644 --- a/src/qsp_codegen/codegen.py +++ b/src/qsp_codegen/codegen.py @@ -61,6 +61,7 @@ def __init__(self, sbml_path: str): self._parse_function_definitions() # before reactions/rules (they may call them) self._parse_reactions() self._parse_rules() + self._parse_rate_rules() # after parameters/species/reactions/rules exist self._parse_initial_assignments() self._parse_events() @@ -231,6 +232,55 @@ def _parse_rules(self): "expression": expr, }) + def _parse_rate_rules(self): + """Promote each on a constant=false parameter to a state. + + SimBiology models an abstract non-physical state (e.g. a dimensionless + niche fraction) as a parameter with ConstantValue=false plus a rate + rule, rather than as a species. The generator has no native rate-rule + path, so we promote the parameter to an amount-tracked pseudo-species + (state slot + SP_ enum + initial value) and inject a single synthetic + production reaction whose rate law IS the rate-rule expression + (d(var)/dt = expr). The existing stoichiometry -> ydot -> Jacobian -> + RHS machinery then evolves it like any other state, with no special + casing downstream. Rate rules whose target is a species (already a + state) are left to the normal reaction path. + """ + params_by_id = {p["id"]: p for p in self.parameters} + for rr in self.model.findall(f".//{SBML_NS}rateRule"): + var_id = rr.get("variable") + param = params_by_id.get(var_id) + if param is None: + continue # target is a species or unknown — not the param path + math = rr.find(f"{MATH_NS}math") + expr = self._mathml_to_infix(math) if math is not None else "0" + var_name = param["name"] + # Promote: drop from parameters, add as a dimensionless amount- + # tracked pseudo-species (has_only_substance_units -> no volume + # scaling; compartment unused). id_to_name[var_id] already maps to + # var_name from _parse_parameters, so reaction-rate references and + # the synthetic product resolve to SPVAR(SP_). + self.parameters = [p for p in self.parameters if p["id"] != var_id] + self.species.append({ + "id": var_id, + "name": var_name, + "base_name": var_name, + "compartment": None, + "initial_value": param["value"], + "is_initial_concentration": False, + "units": param.get("units", "dimensionless"), + "has_only_substance_units": True, + "is_rate_ruled": True, + }) + self.reactions.append({ + "name": f"raterule_{var_name}", + "reactant_ids": [], + "product_ids": [var_id], + "reactant_names": [], + "product_names": [var_name], + "rate_law": expr, + }) + # Supported comparison ops → (C++ op, root convention). # For `lt(a, b)` (a < b), the root function `b - a` crosses 0 upward # when a drops past b → rootsFound=+1 (matches "trigger went false→true"). diff --git a/tests/test_rate_rules.py b/tests/test_rate_rules.py new file mode 100644 index 0000000..4a63f7f --- /dev/null +++ b/tests/test_rate_rules.py @@ -0,0 +1,92 @@ +"""Tests for SBML on a constant=false parameter. + +SimBiology models an abstract non-physical state (e.g. a dimensionless niche +fraction) as a parameter with ConstantValue=false plus a rate rule, rather than +as a species. The generator has no native rate-rule path, so it promotes the +parameter to an amount-tracked pseudo-species and injects a single synthetic +production reaction whose rate law IS the rate-rule expression. The existing +stoichiometry -> ydot -> Jacobian -> RHS machinery then evolves it like any +other state. + +Regression: before this support, such a rate rule was silently dropped — the +parameter stayed frozen at its initial value, with no derivative emitted. +""" +import xml.etree.ElementTree as ET + +from qsp_codegen.codegen import SBMLModel, SBML_NS, MATH_NS + +SNS = SBML_NS.strip("{}") +MNS = MATH_NS.strip("{}") + + +def _build(params_xml: str, rules_xml: str) -> SBMLModel: + doc = ( + f'' + f'' + f"{params_xml}" + f"{rules_xml}" + f"" + ) + root = ET.fromstring(doc) + m = object.__new__(SBMLModel) + m.model = root.find(f"{SBML_NS}model") + m.id_to_name = {} + m.name_to_id = {} + m.function_defs = {} + m.species = [] + m.parameters = [] + m.reactions = [] + m._parse_parameters() + m._parse_rate_rules() + return m + + +_RATE_RULE = ( + '' + f'' + "kf" + "1TLA" + "kdTLA" + "" +) +_PARAMS = ( + '' + '' + '' +) + + +def test_rate_ruled_parameter_promoted_to_state(): + m = _build(_PARAMS, _RATE_RULE) + + # The rate-ruled parameter is removed from parameters (it's now a state); + # the rate constants remain parameters. + assert {p["name"] for p in m.parameters} == {"kf", "kd"} + + # It is added as an amount-tracked dimensionless pseudo-species. + tla = [s for s in m.species if s["name"] == "TLA"] + assert len(tla) == 1 + assert tla[0]["is_rate_ruled"] is True + assert tla[0]["has_only_substance_units"] is True # no volume scaling + assert tla[0]["initial_value"] == 0.0 + + # A single synthetic production reaction (∅ -> TLA) carries the rate-rule + # expression as its rate law, so ydot[TLA] = +rate = the ODE. + rxn = [r for r in m.reactions if r["product_names"] == ["TLA"]] + assert len(rxn) == 1 + assert rxn[0]["reactant_names"] == [] + expr = rxn[0]["rate_law"].replace(" ", "") + assert "kf" in expr and "kd" in expr and "TLA" in expr + + +def test_rate_rule_on_non_parameter_is_left_alone(): + # A rate rule whose target is not a parameter (e.g. a species, already a + # state) is not promoted — the normal reaction path handles it. + rr = ( + '' + f'kf' + ) + m = _build('', rr) + assert m.species == [] + assert m.reactions == [] + assert {p["name"] for p in m.parameters} == {"kf"}