Skip to content

quickstart.md

Latest commit

 

History

History
211 lines (167 loc) · 7.75 KB

quickstart.md

File metadata and controls

211 lines (167 loc) · 7.75 KB
Raw

Quickstart

Load a model and simulate

The block below runs as written after pip install bngsim. It first writes the model it loads, a two-species network in BioNetGen's .net format, so no file is needed beforehand; with a .net of your own, start at from_net.

from pathlib import Path

import bngsim

# A -> B at rate k, in BioNetGen's network (.net) format.
Path("model.net").write_text(
    """\
begin parameters
    1 k 0.5
end parameters
begin species
    1 A() 100
    2 B() 0
end species
begin reactions
    1 1 2 k
end reactions
begin groups
    1 A 1
    2 B 2
end groups
""",
    encoding="utf-8",
)

# Load a BioNetGen .net file
model = bngsim.Model.from_net("model.net")

# Create an ODE simulator
sim = bngsim.Simulator(model, method="ode")

# Run simulation
result = sim.run(t_span=(0, 10), n_points=101)

# Access results as NumPy arrays
print(result.time.shape)          # (101,)
print(result.species.shape)       # (101, 2): one column per species, in model.species_names order
print(result.observables.shape)   # (101, 2): one column per observable, "A" and "B"
print(round(float(result.observables["A"][-1]), 4))  # 0.6738, which is 100 * exp(-5)

Other model formats

.net is the native format. SBML, Antimony and BNGL load through the same Model interface:

sbml_model = bngsim.Model.from_sbml("model.xml")
ant_model  = bngsim.Model.from_antimony("model.ant")  # pip install bngsim[antimony]
bngl_model = bngsim.Model.from_bngl("model.bngl")   # pip install bngsim[bngl]

from_bngl is not a parser. BNGsim simulates reaction networks and BNGL describes rules, so the model block is handed to BNG2.pl generate_network and the resulting .net is loaded through from_net. That makes BNG2.pl a requirement: the bngl extra supplies it, or point $BNGPATH at your own install. Generated networks are cached, so a repeated load of an unchanged model does not re-run BNG2.pl.

See Loading models for the full picture.

Parameter sweeps

# Update parameters in-memory (no file I/O). The next run goes on from the
# state the last one ended in, not from the initial conditions.
model.set_param("k", 1.0)
result2 = sim.run(t_span=(0, 100), n_points=101)

# Reset species to initial conditions
model.reset()

# Set several parameters at once from a dict (this model has one, k)
model.set_params({"k": 0.25})

Stochastic simulation (SSA)

ssa = bngsim.Simulator(model, method="ssa")
result = ssa.run(t_span=(0, 100), n_points=101, seed=42)

Stochastic simulation (PSA)

PSA (partial-scaling approximation) accelerates large-population runs by leaping reaction firings whose minimum reactant population exceeds poplevel = N_c, while preserving Gillespie statistics. The dispatch and the validate_for_ssa gate are shared with SSA, so the same call shape works for .net, SBML, and Antimony models:

import bngsim

# .net
net_model = bngsim.Model.from_net("psa/tcr_signaling.net")
psa = bngsim.Simulator(net_model, method="psa", poplevel=100.0)
result = psa.run(t_span=(0, 300), n_points=301, seed=42)

# SBML / Antimony — same dispatch
sbml_model = bngsim.Model.from_sbml("ssa_sbml/tcr_signaling.xml")
psa = bngsim.Simulator(sbml_model, method="psa", poplevel=100.0)
result = psa.run(t_span=(0, 300), n_points=301, seed=42)

Choose N_c per Lin, Feng, Hlavacek (2019); accuracy degrades when N_c is too small relative to per-species populations. SBML/Antimony models with kineticLaw shapes that SSA cannot handle (e.g. reversible_non_mass_action, non_integer_stoichiometry) raise SsaValidationError at construct time.

Seed semantics for stochastic methods

The seed= keyword on every stochastic entry point follows one contract:

  • seed=None (the default, or omitting seed=): bngsim draws a fresh 31-bit seed from system entropy on each call. Two consecutive calls produce independent trajectories.
  • seed=N (an integer from 0 to 2**64 - 1): bngsim passes N down to the backend verbatim. A fresh Simulator (or NfsimSession / RuleMonkeySession) initialized with the same seed=N reproduces the same trajectory. A seed outside that range raises ValueError.
  • NFsim (method="nf", NfsimSession) seeds its generators from 32 bits, so it takes a seed from 0 to 2**32 - 1 and raises ValueError on a larger one (issue #1515). It used to run such a seed as its low 32 bits, seed=2**32 + 12345 repeating seed=12345, while Result.seed reported the value passed. SSA, PSA and RuleMonkey use all 64 bits.

Result.seed exposes the integer that was actually used, so caller-generated and bngsim-drawn seeds are equally inspectable:

import bngsim

sim = bngsim.Simulator(model, method="ssa")

# Default — fresh seed each call:
r1 = sim.run(t_span=(0, 100), n_points=101)
r2 = sim.run(t_span=(0, 100), n_points=101)
assert r1.seed != r2.seed                        # seeds differ

# Explicit, reproducible:
sim_a = bngsim.Simulator(model.clone(), method="ssa")
sim_b = bngsim.Simulator(model.clone(), method="ssa")
ra = sim_a.run(t_span=(0, 100), n_points=101, seed=42)
rb = sim_b.run(t_span=(0, 100), n_points=101, seed=42)
assert ra.seed == rb.seed == 42                  # exact reproduction

The same contract applies to Simulator.run_batch, NfsimSession.initialize, and RuleMonkeySession.initialize. For run_batch, run_replicates and parameter_scan, the resolved value is the base seed. Row 0 runs at the base itself, so seed=rows[0].seed reruns the whole call, a drawn base included. Rows 1, 2, … take, in order, the distinct values other than the base among the 32-bit words of numpy's SeedSequence(base_seed).generate_state(...), each masked to 31 bits. Every row's seed is stamped onto its Result.seed, so run(seed=r.seed) repeats that row. Bases next to each other give unrelated seeds, and the seeds of one call are distinct. They were base_seed + i, which made an ensemble from seed 1 the one from seed 0 shifted by a row (issue #1365).

Result.seed is None for ODE results and round-trips through HDF5 Result.save() / Result.load(). For sessions, the seed used by initialize() is exposed as session.seed and stamped onto every Result returned by simulate().

Reproducibility unit. The reproducibility contract is:

same starting model state + same seed=N → same trajectory

The Simulator object itself is not the unit. The C++ SSA/PSA backends construct a fresh std::mt19937_64 from the passed seed on every .run() call, so passing the same explicit seed= always seeds the RNG identically. What persists across .run() calls on the same Simulator is the model state (species concentrations) — that's what makes simulate(...); simulate({continue=>1, ...}) multi-segment protocols work. So:

sim = bngsim.Simulator(model, method="ssa")

# First run: starts from model's initial state, seed=42.
r1 = sim.run(t_span=(0, 100), n_points=101, seed=42)

# Second run on the same sim with seed=42 *continues from r1's end
# state*, then advances with the seed=42 RNG stream. r2 is correctly a
# different trajectory because it starts from a different state — not
# because of any RNG issue.
r2 = sim.run(t_span=(0, 100), n_points=101, seed=42)
assert r2.species[0].tolist() == r1.species[-1].tolist()  # state continues

# To re-run the same trajectory: reset the model first.
model.reset()
r3 = sim.run(t_span=(0, 100), n_points=101, seed=42)
assert (r3.species == r1.species).all()  # exactly reproduces r1

The state lives in the Model, not the Simulator: a second Simulator built on the same model also starts from where the last run ended. To run two trajectories from the model's initial state, call model.reset() between them, or give each its own model.clone().