ODE export: design¶

Status: approved in conversation on 2026-10-05, replaces sbmlodefac.

Closes #481 (an updated code generator with complete coverage, run against the test suite) and #438 (initial assignments not supported).

Goal¶

An SBML model is exported as its system of ordinary differential equations in six formats which share one analysis and one math printer core:

• numerical formats, which are code that simulates the model: python (numpy, a scipy OdeSolver), julia (OrdinaryDiffEq.jl of DifferentialEquations.jl), R (deSolve);
• presentation formats, which are documents that describe the model: typst, LaTeX, markdown.

The numerical code is correct, meaning it reproduces roadrunner over the SBML semantic test suite, and readable, meaning a person can read the equations in it. The presentation formats render the same content in the same order, so a model reads the same in every one of them. This is first class functionality of sbmlutils with its own documentation page.

Decisions¶

Why not sbmlmath¶

Issue #481 suggests sbmlmath (MathML to sympy) with the sympy code printers. A probe with sbmlmath 0.4.1 and sympy 1.14.0 showed:

• parsing fails for a piecewise whose condition holds &&;
• rem becomes sympy Mod, which is the floored modulo, while SBML rem has the sign of the dividend;
• log(2, A) is printed by the numpy printer as numpy.log(A, 2), whose second argument is the out array, and by the LaTeX printer as \log(A), dropping the base;
• sympy has no typst printer, and delay has no printer at all;
• sympy evaluates and reorders the expression, so the printed equation is no longer the equation of the model.

Each of these needs an override, so the sympy route ends in custom printers on top of two extra dependencies. A printer on the libsbml AST with a dialect table per language keeps the exact SBML semantics, preserves the structure of the math as written, covers typst, and generalizes python_math, which is already verified against roadrunner.

Architecture¶

sbmlutils/converters/ode/
    __init__.py        public API: OdeSystem, FORMATS, render, write
    system.py          OdeSystem: the frozen dataclass tree of the analysis, from_sbml
    analysis.py        analysis: SBML document to OdeSystem, one method per element
    dependencies.py    ordering of assignments, cycle detection
    events.py          the continuous root function of an event trigger
    astutil.py         shared helpers to build and read the libsbml AST
    printers/
        base.py        MathPrinter: precedence, associativity, dispatch on AST type
        python.py      PythonPrinter
        julia.py       JuliaPrinter
        r.py           RPrinter
        document.py    DocumentPrinter: the base of the typeset dialects
        latex.py       LatexPrinter (markdown uses it inside $$)
        typst.py       TypstPrinter
    symbols.py         naming: code identifiers and typeset symbols
    text.py            free text: single line for code, escaping per document format
    formats.py         Format registry, jinja2 environment, filters, render, write
    documents.py       the context of the document formats, built with the printers
sbmlutils/resources/converters/ode/
    python.py.jinja    julia.jl.jinja    r.R.jinja
    typst.typ.jinja    latex.tex.jinja   markdown.md.jinja

The three layers are independent: the analysis knows no format, a printer knows no model, a template knows no libsbml.

Layer 1: analysis (system.py, analysis.py)¶

OdeSystem.from_sbml(source) takes a path, an SBML string or an SBMLDocument (read with read_sbml) and returns a frozen dataclass tree:

• ModelInfo: id, name, SBML level and version, notes as plain text, model units (time, substance, extent, volume, area, length) as strings from udef_to_string.
• Symbol: id, name, unit string, SBO term, kind (compartment, species, parameter, reaction, species_reference).
• Compartment, Species, Parameter: the symbol plus value, constant, and for a species its compartment, hasOnlySubstanceUnits, boundaryCondition, conversion factor.
• FunctionDefinition: id, arguments, body AST.
• InitialAssignment, Assignment (assignment rule): variable, AST.
• Reaction: id, reactants, products, modifiers with their stoichiometry (a number, or the AST of a stoichiometry given by a rule or an initial assignment of a species reference), reversible flag, rate AST, local parameters. Local parameters are renamed to <reaction id>_<local id> (made unique against all ids) and become constant parameters; the rate AST is rewritten accordingly.
• Ode: state variable, right hand side AST, origin (reactions or rate_rule), and for species in concentration the volume term.
• Event: id, trigger AST, initialValue, persistent, delay AST, priority AST, useValuesFromTriggerTime, assignments (variable, value AST, scale AST, divisor id): the new value is value * scale / new(divisor), where the value is evaluated at trigger or execution time as useValuesFromTriggerTime says, while scale (the size conversion of a species in concentration or held as amount) is always evaluated at execution, and divisor is a compartment assigned by the same event, see the docstring of EventAssignment.
• unsupported: list of (construct, element id), e.g. ("algebraic rule", "rule3").

The following SBML semantics are resolved once, here:

• State variables are the species which are neither constant nor boundary nor assigned and take part in a reaction or a rate rule, plus every parameter, compartment and species reference with a rate rule. A boundary species with a rate rule is a state as well. Every other species, parameter and compartment is constant or assigned.
• Species are stored as roadrunner holds them: in amount if hasOnlySubstanceUnits, else in concentration. The reaction terms of the ODE of a species in concentration are divided by its compartment. A species in concentration whose compartment is not constant (a rate rule, an assignment rule or an event assignment changes its size) is integrated as its amount n_S instead, with the concentration S = n_S / V as an assignment: the amount is what SBML conserves when the size changes, so a rate rule, an assignment rule and an event on the compartment are all exact without a dV/dt term.
• Conversion factors: a species conversion factor, else the model conversion factor, multiplies the reaction terms of the species.
• Stoichiometry: a number, or the id of a species reference with an assignment rule, rate rule or initial assignment, which is then a symbol of the system.
• Rate rules on species in concentration apply to the concentration as written (SBML semantics), no volume term.
• rateOf(x) is replaced by the right hand side of x for a state, by 0 for a constant, and is unsupported for an assigned variable whose derivative would need symbolic differentiation.
• Assignments (assignment rules and the reaction rates) are ordered by their dependencies with graphlib.TopologicalSorter; a cycle raises ValueError naming the variables. Ties keep the order of the document.
• Initial values at t=0: the order of evaluation is the topological order over initial assignments, assignment rules and the initial values of the elements together, so that an initial assignment which depends on an assignment rule (and the reverse) is evaluated right. Initial amounts and concentrations are converted to the representation of the state with the initial size of the compartment, which itself may come from an initial assignment.
• Function definitions are kept as functions, called by the math.
• Unsupported: algebraic rules, delay(), fast reactions, an event assignment to a constant, the comp, fbc and distrib packages (a comp model is flattened with flatten_sbml first). They are collected, never silently dropped.

Layer 2: math printers (printers/)¶

MathPrinter.print(ast, symbols) returns a string. The base class implements the traversal: dispatch on ASTNode.getType(), parentheses from a precedence and associativity table, n-ary operators, unary minus, and the structure of piecewise, log with base, root with degree, relations with more than two operands. A dialect overrides tables and a few node handlers:

Booleans are numbers in SBML: a relation used as a number becomes float(...) in python, Float64(...) in julia, as.numeric(...) in R. The dialects follow IEEE 754 for a relation with NaN (false, and true for !=): R, where such a relation is NA, writes every relation as isTRUE(a > b) and a != b as !isTRUE(a == b). Julia writes the functions whose Base version throws outside its real domain (sqrt, log, asin, acosh, a power with a negative base and a fractional exponent, ...) with NaNMath, so the result is NaN as in SBML; every julia number is a Float64 literal, so that integer arithmetic cannot overflow. Every construct the dialect cannot express raises NotImplementedError naming the construct. Identifiers come only from the symbols mapping; an identifier without a mapping raises, so an id is never written raw.

Presentation printers break long equations: the top level sum of a right hand side is split into lines of at most width terms (default 4) inside an align/aligned block.

Layer 3: formats (formats.py, templates)¶

FORMATS maps a name to a Format: template, printer, file suffixes, kind (code or document) and the options it accepts. The API:

from sbmlode import OdeSystem

system = OdeSystem.from_sbml("model.xml")
code: str = system.render("python", simulator=True)
system.write("model.py")  # format from the suffix
system.write("model.typ", standalone=False)
system.render_template(Path("my_template.jinja"))

render checks the options against the format, builds the symbol mapping of the format (symbols.py), prints every AST once with the printer of the format and renders the template with a context of plain strings and lists, so a template never touches libsbml. A numerical format raises NotImplementedError listing system.unsupported if it is not empty; a document lists them in a section of their own.

Safety is kept from odefac: every free text goes through single_line in code comments and through the escaping of the format in documents; every id written into code is checked to be an SId; generated python names avoid keywords and the names the code uses itself (t, x, p, np, ...), likewise for julia and R.

Numerical formats¶

The three files have the same structure; python is 0-indexed, julia and R are 1-indexed, and the indices are hidden behind named locals.

1. Header: model id and name, source, sbmlutils version, units of the model, the unsupported constructs (empty for a file that is written at all).
2. Id tables xids, pids, yids, one entry per line with name and unit as a comment.
3. p: default values of the constant parameters, constant compartments and constant species.
4. initial_values(p): returns (x0, p), the constants set by an initial assignment updated in p, evaluating initial values, initial assignments and assignment rules at t=0 in their order. This fixes #438.
5. Right hand side: f_dxdt(t, x, p) in python (the order of the scipy solvers), f!(dx, x, p, t) in julia (every julia function takes (x, p, t), the order of DifferentialEquations.jl), f_dxdt(t, x, p) returning list(dx) in R. The body unpacks named locals from x and p, calls the function definitions (emitted as functions before it), evaluates the assignments in order, the reaction rates under their reaction ids, then one line per state, each with the name and unit as a comment. f_y(t, x, p) returns all assigned values and reaction rates.
6. Events (simulator=True and False): event_triggers(t, x, p) returns one continuous root function value per event and event_conditions(t, x, p) the exact truth value of every trigger; event_values_<id>(t, x, p) evaluates the assigned values (at trigger or execution time as useValuesFromTriggerTime says) and event_assign_<id>(t, x, p, values) applies them at execution, with the size conversions scale/divisor evaluated then; EVENTS lists per event its flags, delay and priority functions. A single relation a > b becomes a - b, && the minimum, || the maximum and ! the negation of its operands' root functions, so that the sign of the root function is the truth value of the trigger. Strict and non-strict relations differ only at the root, which the event handling resolves with event_conditions.
7. simulate(t_end, points=101, p=None, x0=None, ...) (only simulator=True; python adds rtol, atol, method, max_step, max_steps): integrates with a limit of solver steps and a guard against a step size which no longer advances the time, both raising, so a model which blows up fails fast, stops where a trigger changes its truth value, evaluates triggers with initialValue, orders simultaneous events by priority, schedules delayed events, drops non-persistent events whose trigger turned false before execution, applies assignments with the values of trigger or execution time, re-evaluates the assigned variables, restarts the integration, raises after a limit of cascaded executions at one time (10000, as roadrunner), and returns a table with time, the states, the assigned values and the constants changed by events (pandas DataFrame, DataFrames.DataFrame, data.frame). Python steps a scipy.integrate OdeSolver (LSODA, rtol=1e-8, atol=1e-10 as defaults, overridable) and finds a change of a trigger by bisection on the dense output with the exact boolean triggers; julia steps an OrdinaryDiffEq integrator (init/step!, Rodas5P) with a VectorContinuousCallback ending a step at a root, R integrates segment by segment with deSolve::lsoda and rootfunc in an own loop (deSolve events know no strict relations, delays, priorities or persistence); julia and R reuse the execution algorithm of python (execute_events) and resolve strict versus non-strict relations at a root with event_conditions. The python file runs as a script (if __name__ == "__main__") and prints the head of the table.

The correctness contract: initial_values, f_dxdt, f_y and simulate match roadrunner (rtol=1e-6, atol=1e-9) for every case of the SBML semantic test suite without an unsupported construct.

Presentation formats¶

One section order for all three formats, the sections without content are omitted:

1. Title (model name, else id), a metadata line (id, SBML level and version, source, sbmlutils version), the notes of the model as plain text.
2. Units of the model.
3. Compartments, species and parameters: tables with symbol, id, name, value, unit and the constant flag; species add compartment, amount or concentration and boundary condition.
4. Function definitions as f(x, y) = ....
5. Initial assignments and assignment rules, in their order.
6. Reactions: the reaction equation (2 A + B -> C, <-> if reversible, modifiers after it), the rate v_r = ... and the local parameters.
7. ODE system: d x / d t = ... per state, written with the reaction rates (v_1 - 2 v_2), conversion factors and volumes explicit, rate rules marked.
8. Events: trigger, delay, priority, flags and assignments.
9. Unsupported constructs.

Symbols: an id is typeset as a math symbol. The part before the first _ is the base, the rest the subscript (k_cat_glc is k with subscript cat_glc); a base of more than one letter is upright. The option symbols="id" | "name" selects ids or element names (a name which is no valid symbol falls back to the id).

Format specifics:

• typst: standalone=True writes #set document, page and text settings, headings and the content; False the content only. Math in $ ... $, tables with #table. Escaped characters in text: backslash, #, $, *, _, @, <, >, [, ] and the backtick.
• LaTeX: standalone=True writes an article with amsmath, amssymb, booktabs, xltabular, parskip and hyperref, the tables being longtable, or xltabular where a column of long ids must wrap, and iftex choosing fontenc and lmodern for pdfLaTeX and fontspec for XeLaTeX and LuaLaTeX, so that it compiles with every engine; False the body only and a comment listing the packages it needs. Escaping by tex_text.
• markdown: GitHub flavored tables, display math in $$ ... $$ with the LaTeX printer (rendered by GitHub and by Zensical with MathJax), ids in code spans. Escaped: |, *, _, backtick, <, > and the HTML entities. This replaces the markdown which create_model(create_markdown=True) writes.

Testing¶

Tests live in tests/converters/ode/, as test_ode_<topic>.py with the shared helpers in ode_helpers.py (the directory has no __init__.py, the unique basenames keep the import mode of pytest working):

• test_ode_system.py: the analysis, one test per semantic rule above (conversion factors, variable compartments, local parameter renaming, stoichiometry math, rate rules on every kind, rateOf, ordering, cycles, initial values, unsupported collection); test_ode_astutil.py, test_ode_symbols.py and test_ode_text.py the helpers of the analysis, the naming and the text escaping.
• test_ode_printers.py: every AST node type in every dialect against golden strings, the precedence table, and the evaluation of the numerical dialects against roadrunner for a set of formulas which covers every construct of the math.
• test_ode_python.py, test_ode_julia.py and test_ode_r.py: compare initial_values, f_dxdt, f_y and the trajectory of simulate with roadrunner. Python always (skips without roadrunner); julia and R are run through the command prefixes in SBMLUTILS_JULIA and SBMLUTILS_RSCRIPT (default julia and Rscript, split like a shell command, so they can be a docker run) and skip when the toolchain or its packages are missing, unless SBMLUTILS_REQUIRE_TOOLCHAINS=1 makes that a failure. The tolerances against roadrunner are rel=1e-8 for the values at a point and rtol=1e-6, atol=1e-9 for a trajectory, relaxed to 1e-4 and 1e-6 for a model with events, whose times are located to the tolerance of the integration.
• test_ode_presentation.py: typst compiled with the typst python package (added to the dev extra), LaTeX compiled with tectonic when it is on the path, markdown parsed with markdown-it-py; golden files of the demo model, the repressilator and a model with events in tests/converters/ode/golden/, regenerated with the environment variable SBMLUTILS_UPDATE_GOLDEN=1.
• test_ode_safety.py: the injection and SId tests, run on every format.
• test_ode_docs.py: the files of docs/images/ode, which docs/ode.md shows, are the current output of the export.
• test_ode_testsuite.py: a curated subset of the SBML semantic test suite (CURATED, every feature tag at least once) runs in the default test run for python and, with the toolchain, in the tox environments for julia and R; the full sweeps run behind the sbml_testsuite marker, a known failure is a strict xfail with its reason. scripts/ode_report.py runs the sweep, one process per case, and reports the pass rate per numerical format with the known failures and their reasons, as scripts/roundtrip_report.py does.

Continuous integration: tox environments julia and R (python 3.14 plus the toolchain), jobs in ci-cd.yml with julia-actions/setup-julia (the packages of tests/converters/ode/julia/Project.toml: OrdinaryDiffEq, DataFrames, NaNMath, SpecialFunctions) and r-lib/actions/setup-r (deSolve, from binary packages), and a latex job with tectonic. They are not required checks of the ruleset at first, so a toolchain outage does not block a merge; the decision to make them required is left to the maintainer.

Documentation¶

• docs/ode.md, "ODE export": the API, the options, one section per format with the output of the repressilator (code listings for the numerical formats, the typst output as SVG pages, the LaTeX source, the markdown rendered), how to run the generated code with its solver, and the table of supported SBML features. The SVG pages of the typst document are not committed: examples/converters/ode.py compiles them (python -m examples.converters.ode docs/images/ode) with the typst package, deterministically with the fonts of typst only, which the documentation workflow runs before the build. The text files it writes next to them are committed, shown through pymdownx.snippets, and kept current by test_ode_docs.py. MathJax, which typesets the math of the markdown output, is vendored in docs/javascripts/mathjax/.
• docs/api/converters.ode.md replaces docs/api/converters.odefac.md; docs/converters.md links to the new page.
• examples/converters/ode.py writes all six formats of a packaged model.

Migration¶

sbmlodefac and its templates are removed. create_model(create_markdown=True), examples/assignment.py and examples/tutorial/linear_chain.py use the new API. The release notes of the next release map the old calls to the new ones: