Conservation laws & parameter scans

Conservation laws

BNGsim automatically detects conservation laws at model load time via Gaussian elimination on the stoichiometry matrix. A conservation law is a linear combination of species that remains constant during the dynamics:

Σ L[k,i] · y[i] = constant_k    for all time

Conservation laws arise from molecular conservation (total receptor, total ligand, etc.) and are detected for ALL input formats (.net, Antimony, SBML, programmatic ModelBuilder).

model = bngsim.Model.from_net("model.net")

# Inspect conservation laws
laws = model.conservation_laws
print(laws["n_laws"])           # number of conservation laws
print(laws["dependent"])        # dependent species indices (0-based)
print(laws["independent"])      # independent species indices
print(laws["constants"])        # conservation constants from ICs
print(laws["coefficients"])     # n_laws × n_species coefficient matrix

Impact on the steady-state solver: Models with conservation laws have a rank-deficient Jacobian, which causes standard Newton solvers to fail. BNGsim’s reduced-space Newton solver automatically handles this:

  1. Identifies N - n_laws independent species from the conservation structure

  2. Builds a reduced residual function on the independent subspace

  3. Solves the non-singular reduced system via KINSOL

  4. Reconstructs dependent species from the conservation constraints

This is transparent to the user — sim.steady_state() just works, regardless of whether the model has conservation laws.

parameter_scan integration

BNGsim supports BioNetGen’s parameter_scan action, commonly used in BNGL models for dose-response analysis. When PyBNF encounters a parameter_scan action in a BNGL model, it automatically routes through BNGsim’s batch steady-state or time-course infrastructure.

For BNGL models that include:

parameter_scan({method=>"ode", parameter=>"L_0", \
    par_min=>0.01, par_max=>100, n_scan_pts=>50, \
    log_scale=>1, steady_state=>1, t_end=>1e6})

PyBNF’s BngsimModel parses this action and dispatches to:

  • strict BNG2.pl parity (default) when steady_state=>1: each scan point runs sim.run(steady_state=True) — the same run_network -c integrate-to-||f||_2/n early-stop BNG2.pl uses.

  • sim.steady_state_batch() (KINSOL Newton, parity fallback) when the model opts in with ss_method=>"newton" (or its alias ss_method=>"kinsol") — a speed accelerator for monostable dose-response sweeps.

  • Time-course simulation to t_end when steady_state=>0.

ss_method=>"newton" is rejected for bifurcate continuation scans (warn + downgrade to the parity path): bifurcate carries state between points to detect hysteresis/multistability, and independent-per-point Newton finds a root — it can jump branches and destroy the hysteresis signal.

The output is a .scan file with columns for the scanned parameter and all observables/expressions, matching BNG2.pl’s parameter_scan output format. This enables transparent acceleration of existing BNGL workflows that use parameter_scan for dose-response curve generation.