simulation.md

Running simulations

Batch parameter sweeps

import numpy as np

# Sweep over parameter values (GIL released, thread-parallel)
param_sets = [{"kf": v} for v in np.linspace(0.1, 10, 100)]
results = sim.run_batch(
    t_span=(0, 100),
    n_points=101,
    params=param_sets,
    num_processors=4,       # parallel threads
)
# results is a list of 100 Result objects

# Or get a single Result with 3D arrays
batch = sim.run_batch(
    t_span=(0, 100), n_points=101,
    params=param_sets, squeeze=True,
)
print(batch.species.shape)  # (100, 101, n_species)

run_batch clones the model and resets every point to the .net seed initial conditions — correct for a dose–response from a fresh start. Where the model is not at them (after a set_concentration, or a run), it warns that the points do not start where the model is. run_replicates starts every replicate where the model is, as run does.

Parameter scans that carry a pre-equilibrated state

A run_batch re-syncs each point to the seed, which is wrong for a multi-phase protocol: pre-equilibrate, intervene (wash / bolus), then scan from the post-intervention state. parameter_scan matches BNG2.pl's reset_conc semantics — each point resets to the state at scan invocation (or to a named snapshot), not the seed:

# Phase 1 — pre-incubate off the seed, then snapshot the post-wash state.
sim.run_until(t=7200)                    # load receptors
sim.model.set_concentration("L(r,label~hot)", 0.0)   # wash out free hot ligand
sim.save_concentrations("start_competition")

# Phase 2 — dose-response that resets each point to the saved snapshot.
results = sim.parameter_scan(
    "cold_conc",
    par_scan_vals=[2.3e-10, 6.6e-9, 3.9e-8, 4.9e-7],
    t_span=(0, 1200), n_points=13,
    reset_conc=True,                     # BNG reset_conc=>1
    reset_to="start_competition",        # reset each point to the snapshot
    # Apply a coupled setConcentration that tracks the scanned value:
    on_point=lambda model, v: model.set_concentration(
        "L(r,label~cold)", v * NA * Vecf
    ),
)
# Each Result carries custom_attrs["scan_value"]; the bound-hot readout starts
# near the pre-wash amount and falls as cold competes it off.

# bifurcate() is the continuation sibling (reset_conc=>0): each point continues
# from the previous point's end-state — sweep up then down to trace hysteresis.
branch = sim.bifurcate("stimulus", par_min=0, par_max=10, n_scan_pts=50,
                       t_span=(0, 500), n_points=2)

Interactive simulation

# Run to a time point, inspect/modify, continue
sim.run_until(t=50)

# Knock out a reaction mid-simulation
sim.intervene({"kf": 0.0})

# Continue from current state
result = sim.run_until(t=100)

# Save and restore snapshots
snap = sim.snapshot()
sim.run_until(t=200)
sim.restore(snap)  # back to t=50

A leg continues the events as well as the species and the clock. Each trigger keeps the truth it had when the previous leg ended, so an event whose trigger is still true at a boundary does not fire again there, and a delayed execution pending across the boundary is applied in the next leg (issue #693). An intervention that makes a trigger true between legs fires it at the start of the next leg. snapshot() and restore() save and rewind the events with the state, and so does rolling the clock back to where one of the last 64 legs of this trajectory ended with set_state(x, time=t) or set_time(t) (below); rolling back further warns and starts the events afresh, so take a snapshot() there instead. A run from any other start time, or after model.reset(), is a fresh start: each trigger starts from its initialValue.

Rolling a step back

restore rebuilds the backend, which is the right cost for an occasional checkpoint and the wrong one inside a loop. A predictor-corrector coupling loop redoes every step: advance, correct the coupling input, advance the same interval again from the saved point. set_state(state, time=t) moves the state and the clock together and keeps the warm ODE solver:

saved, t0 = sim.get_state().copy(), sim.current_time
sim.run_until(t0 + dt)                   # predictor
corrected = saved.copy()
corrected[idx] = coupling_value          # e.g. from the other half of the split
sim.set_state(corrected, time=t0)        # back to t0, with the correction
sim.run_until(t0 + dt)                   # corrector, over the same interval

The clock matters whenever the model reads it — a time-indexed table function, time() in a rate law or function. Rolling back the state alone leaves current_time a step ahead, and the redone step then sees its forcing at the wrong times. sim.set_time(t) sets the clock on its own. Both calls validate before they assign: a non-finite time raises ValueError, and a rejected call changes neither the state nor the clock.

ReactionKernel.set_state(state, time=t) does the same for a kernel.

A run_until step continues the trajectory the previous one ended, which the ODE solver uses to keep its linear-solver state across the steps of a loop (see Linear solver). An independent run() starts afresh.

Species manipulation

# Save current concentrations as new initial state
model.save_concentrations()

# Perturb a single species (e.g., bolus addition)
model.set_concentration("A(b)", 500.0)

# Check concentration
print(model.get_concentration("A(b)"))

# Reset restores to save_concentrations() snapshot
model.reset()

# Named states coexist — save two, restore either (BNG saveConcentrations("name"))
model.save_concentrations("t=0")
# ... advance / intervene ...
model.save_concentrations("start_competition")
model.restore_concentrations("t=0")            # back to the first named state
model.restore_concentrations("start_competition")

Stop conditions

# Expression-based: stop when observable drops below threshold
sim.add_stop_condition("A_total < 10", label="low_A")

# Callable-based: arbitrary Python logic
sim.add_stop_condition(
    lambda r: r.species[-1, 0] < 5.0,
    label="very_low_A",
)

try:
    result = sim.run(t_span=(0, 1000), n_points=1001)
except bngsim.StopConditionMet as e:
    print(f"Stopped at t={e.result.time[-1]}: {e.condition}")
    partial = e.result  # truncated result up to stop point

A condition is checked against the completed result, so it does not interrupt the solve. The backend integrates the whole interval and leaves the model at its end; e.result is that trajectory truncated at the first row where the condition held, and sim.current_time reads the end of the interval rather than the trigger time. The clock is the one that names the state the model holds, so a following run_until continues from there correctly.

To continue from the trigger point instead, snapshot first:

snap = sim.snapshot()
try:
    sim.run_until(t=1000)
except bngsim.StopConditionMet as e:
    sim.restore(snap)   # back to where the leg started

run_batch, run_replicates and parameter_scan (with bifurcate) check every row the same way and raise nothing: a row that meets a condition comes back cut there, beside the rows that ran to the end, and Result.stop_condition names the condition it met (None on a row that met none). Conditions are checked in the order they were added, and the first one met cuts the row. squeeze=True is refused while a condition is registered, since a cut row is shorter than the rest; Result.squeeze(rows) stacks the list when no row was cut.

rows = sim.run_batch(params=[{"k": 0.1}, {"k": 0.001}], t_span=(0, 100))
[r.stop_condition for r in rows]   # e.g. ['low_A', None]