Stochastic simulation#

Given a grammar, Stochsim runs an exact Gillespie (stochastic simulation algorithm, SSA) simulation. Rather than enumerating the whole reaction network up front, it grows it on demand as new species appear, while evolving the molecule counts one reaction event at a time. You supply the initial counts and a per-reaction rate; simulate() advances the trajectory.

Two idealized reactors bracket the behavior. A continuous stirred-batch reactor (CSBR) is closed — no material flows in or out — and relaxes to a steady state; a continuous-flow stirred-tank reactor (CFSTR) is open, with in- and out-flows, and can reach non-equilibrium steady states. The example below is a closed batch.

Atrazine degradation#

The anaerobic soil degradation of the herbicide atrazine to cyanuric acid is driven by two reaction types — hydrolysis and reductive dealkylation. The triazine ring is modeled with one labeled node per substituent site, so each rule is a single site relabeling written in DFS notation (ruleDFS("[C]1[Cl]2>>[C]1[OH]2"), and likewise Et → Am, iPr → Am, Am → OH). In a closed batch every intermediate can still react, so the system deadlocks only once every molecule has become cyanuric acid:

def triazine(s0, s2, s4, name):
    return graphGMLString(f'''graph [
      node [ id 0 label "C" ] node [ id 1 label "N" ] node [ id 2 label "C" ]
      node [ id 3 label "N" ] node [ id 4 label "C" ] node [ id 5 label "N" ]
      edge [ source 0 target 1 label "-" ] edge [ source 1 target 2 label "-" ]
      edge [ source 2 target 3 label "-" ] edge [ source 3 target 4 label "-" ]
      edge [ source 4 target 5 label "-" ] edge [ source 5 target 0 label "-" ]
      node [ id 6 label "{s0}" ] edge [ source 0 target 6 label "-" ]
      node [ id 7 label "{s2}" ] edge [ source 2 target 7 label "-" ]
      node [ id 8 label "{s4}" ] edge [ source 4 target 8 label "-" ] ]''', name, add=False)

rules = [ruleDFS(f"[C]1[{a}]2>>[C]1[{b}]2", rid, add=False) for a, b, rid in [
    ("Cl", "OH", "halogen_hydrolysis"), ("Et", "Am", "de_ethylation"),
    ("iPr", "Am", "de_isopropylation"), ("Am", "OH", "amine_hydrolysis")]]

def rate(he):
    n = {r.name for r in he.rules}
    return 3.32e-8 if "de_ethylation" in n else 2.65e-8 if "de_isopropylation" in n else 5.00e-9

atrazine = triazine("Cl", "Et", "iPr", "atrazine")
cyanuric = triazine("OH", "OH", "OH", "cyanuric")
sim = Stochsim(graphDatabase=[atrazine], expandStrategy=rules,
               initialState={atrazine: 100}, reactionRate=rate, seed=2024)
sim.simulate(iterations=100_000)
print(f"after {sim.iteration} events (t = {sim.time:.3g}): "
      f"atrazine {sim.state(atrazine)}, cyanuric acid {sim.state(cyanuric)} of 100")

after 500 events (t = 1.47e+09): atrazine 0, cyanuric acid 100 of 100

Read the trajectory back as data with state(), series() and trajectory() (times plus per-species counts). A rate callback may return a constant or a (rate, cache) pair; unless it caches, it is re-evaluated each step, so rates can depend on the current time or state(...). print() writes a standalone, seekable interactive HTML playback of the run to the summary folder.

17_stochsim_atrazine.py#

"""17. Rule-based stochastic simulation: atrazine degradation.

Given a grammar, ``Stochsim`` runs an exact Gillespie (SSA) simulation, growing
the reaction network on demand as new species appear. You supply the initial
molecule counts and a per-reaction rate; ``simulate()`` advances the trajectory.

The chemistry is the anaerobic soil degradation of the herbicide atrazine to
cyanuric acid, driven by two reaction types (hydrolysis and reductive
dealkylation). The triazine ring is modeled with one labeled node per
substituent site, so a rule is a single site relabeling:

    halogen hydrolysis:  Cl  -> OH
    de-ethylation:       Et  -> Am   (ethylamino   -> amino)
    de-isopropylation:   iPr -> Am   (isopropylamino -> amino)
    amine hydrolysis:    Am  -> OH

Atrazine {Cl, Et, iPr} drains to cyanuric acid {OH, OH, OH}. In a closed system
every intermediate can react further, so the system deadlocks only when every
molecule has become cyanuric acid.

Run:  python 17_stochsim_atrazine.py
"""
import _bootstrap
from molfoundry import (
    graphGMLString, ruleDFS, DG, Stochsim, GraphPrinter,
    addSubset, repeat, enable_post,
)

enable_post()

# Rate constants (Erickson et al.); hydrolysis shared by halogen & amine.
K_HYDROLYSIS = 5.00e-9
K_DEETHYLATION = 3.32e-8
K_DEISOPROPYLATION = 2.65e-8


def triazine(s0, s2, s4, name):
    """A triazine ring C-N-C-N-C-N with three substituent-site labels."""
    return graphGMLString(
        f'''graph [
          node [ id 0 label "C" ] node [ id 1 label "N" ] node [ id 2 label "C" ]
          node [ id 3 label "N" ] node [ id 4 label "C" ] node [ id 5 label "N" ]
          edge [ source 0 target 1 label "-" ] edge [ source 1 target 2 label "-" ]
          edge [ source 2 target 3 label "-" ] edge [ source 3 target 4 label "-" ]
          edge [ source 4 target 5 label "-" ] edge [ source 5 target 0 label "-" ]
          node [ id 6 label "{s0}" ] edge [ source 0 target 6 label "-" ]
          node [ id 7 label "{s2}" ] edge [ source 2 target 7 label "-" ]
          node [ id 8 label "{s4}" ] edge [ source 4 target 8 label "-" ]
        ]''',
        name,
        add=False,
    )


def site_rule(a, b, rid):
    """Relabel one substituent site a -> b on a ring carbon.

    Written in DFS line notation: "[C]1[a]2>>[C]1[b]2" keeps the ring carbon
    (id 1) and relabels its substituent (id 2) from a to b. The multi-character
    site labels are bracketed, so they are opaque labels rather than atoms.
    """
    return ruleDFS(f"[C]1[{a}]2>>[C]1[{b}]2", rid, add=False)


rules = [
    site_rule("Cl", "OH", "halogen_hydrolysis"),
    site_rule("Et", "Am", "de_ethylation"),
    site_rule("iPr", "Am", "de_isopropylation"),
    site_rule("Am", "OH", "amine_hydrolysis"),
]


def reaction_rate(he):
    names = {r.name for r in he.rules}
    if "de_ethylation" in names:
        return K_DEETHYLATION
    if "de_isopropylation" in names:
        return K_DEISOPROPYLATION
    return K_HYDROLYSIS  # halogen or amine hydrolysis


atrazine = triazine("Cl", "Et", "iPr", "atrazine")
cyanuric = triazine("OH", "OH", "OH", "cyanuric-acid")

# 1. The full network (built eagerly, just to report its size).
full = DG(graphDatabase=[atrazine])
full.build().execute(addSubset([atrazine]) >> repeat(rules), verbosity=0)
print(f"full degradation network: {full.numVertices} species, {full.numEdges} reactions")

# 2. Stochastic simulation of a closed batch of 200 atrazine molecules.
N = 200
sim = Stochsim(
    graphDatabase=[atrazine],
    expandStrategy=rules,
    initialState={atrazine: N},
    reactionRate=reaction_rate,
    seed=2024,
)
sim.simulate(iterations=100_000)

print(f"\nsimulated {sim.iteration} reaction events, clock at t = {sim.time:.3g}")
print(f"atrazine remaining: {sim.state(atrazine)}   "
      f"cyanuric acid formed: {sim.state(cyanuric)}   (started with {N})")

# 3. The time course as data: sample the atrazine and cyanuric-acid counts.
traj = sim.trajectory()
atr_series = sim.series(atrazine)
cyan_series = sim.series(cyanuric)
print("\n   t            atrazine  cyanuric")
n_steps = len(traj.times)
for i in range(0, n_steps, max(1, n_steps // 8)):
    print(f"   {traj.times[i]:.3e}   {atr_series[i]:6d}   {cyan_series[i]:6d}")

# 4. An interactive, seekable HTML playback of the run (a MolFoundry extra).
sim.print(GraphPrinter(), name="atrazine_run")
print("\nWrote summary/atrazine_run.html (interactive trajectory playback).")