chemart / srsim

hub:
  repo_type: generator
  builtin: srsim
  tags:
  - banzhaf-yamamoto
chemistries:
- id: srsim
  name: SRSim (rule-based spatial simulator)
  intuition: 'Rule-based chemistry that refuses to forget geometry. Molecules are spheres carrying binding
    sites at fixed positions, diffusing in three dimensions, and rules are written over patterns as in
    any rule-based language - but a bimolecular rule fires only when the two sites are actually close
    enough and correctly oriented. Geometry therefore decides which complexes can exist at all: change
    a monomer''s binding angle and you get rods instead of rings. The species set is never enumerated,
    only explored.

    '
  aliases:
  - spatial rule-based modelling
  - BNGL + LAMMPS
  origin: Gruenert & Dittrich, 2010-2013
  book: 18.3.3
  refs:
  - '[351]'
  - '[388]'
  - '[415]'
  - '[861]'
  - doi:10.1186/1471-2105-11-307
  family: systems-biology
  kind: framework
  constructive: true
  S:
    repr: 'elementary molecules (EM): spheres with a mass, a radius and a diffusion coefficient, carrying
      named components (binding sites) at fixed polar coordinates (dist, theta, phi) and optional internal
      states. A species is a complex: a connected site graph, taken up to isomorphism, whose id and structure
      are its canonical text, e.g. M(a[.],b[1]),M(a[1],b[.]). Because bonds are the straight connection
      of two component vectors, a complex also has a shape.

      '
    definition: implicit
    notes: 'The species set is open: complexes are created on the fly by the rules, never enumerated in
      advance, which is the book''s point about combinatorial explosion. Complexes whose canonical text
      exceeds 120 characters get a formula plus a hash of that text as id.

      '
  R:
    definition: implicit
    arity:
    - 1
    - 2
    scheme: 'BNGL-style rules over molecule patterns, e.g. ''grow'' M(a) + M(b) <-> M(a!1).M(b!1) @ k,
      k_rev. A rule forms one bond, breaks one bond or changes one internal state. A bimolecular rule
      fires only between molecules that are geometrically compatible: their distance is within distance_tolerance
      of the ideal bond length d_ij + d_kl and, for every bond either molecule already carries, the angle
      to the partner is within angle_tolerance of the ideal angle computed from the component coordinates.
      A compatible pair then reacts with probability 1 - exp(-k dt).

      '
    notes: 'Geometry, not the rule set, decides which complexes exist: the same polymerisation rule gives
      rods, rings or helices depending on the angle between the monomer''s two components (paper fig.
      3), and complexes the geometry forbids are never built. Reactions inside one complex are treated
      as bimolecular, as in the paper''s fig. 1d.

      '
  A:
    reactor: continuous-space
    dilution: none
    notes: 'Coarse-grained molecular dynamics in a periodic cuboid box: harmonic bonds and angles, a soft-sphere
      repulsion and Langevin (or overdamped Brownian) motion, with the rule system evaluated after every
      position update. A broken bond makes both molecules refractory for refractory_time, which suppresses
      geminate recombination. well_mixed=true replaces the space by the well-stirred limit (a pair is
      compatible with probability V_react/V_reactor), which is the paper''s non-spatial BioNetGen control.

      '
  params:
  - name: model
    type: enum
    default: polymer
    choices:
    - polymer
    - dimerization
    - scaffold
    - custom
    role: structural
    meaning: 'polymer: one monomer with two binding sites whose ideal angle is bond_angle (the paper''s
      fig. 3 study of geometry-determined assembly); dimerization: A + B <-> AB, the reversible reaction
      the paper''s supplement analyses; scaffold: the paper''s scaffold-protein model (figs. 5-6) with
      the geometry of its spass.geo and the rates of its spass_w.bngl; custom: molecule_types, rules and
      init'
  - name: n_molecules
    type: int
    default: 50
    min: 1
    max: 2000
    range: the paper's scaffold run uses 300 A and 30 S in a 100^3 box for 100000 steps
    role: population
    meaning: 'copies of the main molecule: M (polymer), A and B each (dimerization), A (scaffold)'
  - name: n_scaffolds
    type: int
    default: 10
    min: 0
    max: 500
    range: the paper uses 30
    role: population
    meaning: 'scaffold model only: copies of the scaffold protein S, which has four binding sites'
  - name: bond_angle
    type: float
    default: 90.0
    min: 0.0
    max: 180.0
    range: 'fig. 3: 180 gives rods, 90 closed quadratic structures, intermediate angles rings and helices'
    role: spatial
    meaning: 'polymer model only: ideal angle between the monomer''s two binding sites, which is what
      decides the shape of the assemblies'
  - name: k_on
    type: float
    default: 0.5
    min: 0.0
    role: kinetic
    meaning: 'polymer and dimerization: microscopic binding rate k2_mic, applied to a geometrically compatible
      pair (probability 1 - exp(-k dt) per step)'
  - name: k_off
    type: float
    default: 0.01
    min: 0.0
    role: kinetic
    meaning: 'polymer and dimerization: microscopic unbinding rate of a bond'
  - name: scaffold_binding
    type: float
    default: 1.0
    min: 0.0
    max: 1000.0
    range: figure 6 compares 1 ('w-Scf', the published d2_on = 0.05) with 0 ('wo-Scf')
    role: kinetic
    meaning: 'scaffold model only: multiplier of the published scaffold association rates d2_on and d3_on,
      i.e. the rate ks of figure 6'
  - name: rate_scale
    type: float
    default: 1.0
    min: 0.0
    role: kinetic
    meaning: multiplies every rule rate, like SRSim's preFactBindR/preFactBreakR/preFactModifyR prefactors;
      scaling the whole rule set rescales time, so rate_scale x simulated time reproduces a much longer
      published run
  - name: box
    type: float
    default: 12.0
    min: 1.0
    range: the published examples use 80^3, 100^3 and 200x30x30
    role: spatial
    meaning: edge of the cubic reaction volume; the box is periodic and must be wider than twice the largest
      reaction distance
  - name: steps
    type: int
    default: 1500
    min: 1
    max: 200000
    range: the published runs are 100000 to 3000000 steps
    role: kinetic
    meaning: molecular-dynamics steps; reactions are looked for after every one of them
  - name: dt
    type: float
    default: 0.02
    min: 1.0e-09
    range: 'the paper''s rule of thumb: sqrt(6 D dt) should be about a tenth of a particle diameter'
    role: kinetic
    meaning: time step of the molecular dynamics
  - name: diffusion
    type: float
    default: 1.0
    min: 0.0
    range: the paper estimates 8e-11 m^2/s for a haemoglobin-sized protein (Stokes-Einstein)
    role: spatial
    meaning: diffusion coefficient of an elementary molecule (the scaffold S gets a fifth of it, being
      five times heavier)
  - name: temperature
    type: float
    default: 1.0
    min: 1.0e-09
    role: thermodynamic
    meaning: k_B T of the Langevin thermostat; the friction is gamma_0 = k_B T / D
  - name: distance_tolerance
    type: float
    default: 0.5
    min: 0.0
    range: every published example uses DistanceDeviation = 0.4 or 0.5
    role: spatial
    meaning: how far the distance of two molecules may deviate from the ideal bond length for them to
      react
  - name: angle_tolerance
    type: float
    default: 30.0
    min: 0.0
    max: 180.0
    range: 'published AngularDeviation: 180 (scaffold), 50 (microtubule and spheres), 30 (Sierpinski tiles)'
    role: spatial
    meaning: how far the angles at an already bound molecule may deviate from their ideal values for a
      new bond to be accepted; 180 accepts any direction
  - name: refractory_time
    type: float
    default: 1.0
    min: 0.0
    range: the published runs use 50 to 200 time steps
    role: kinetic
    meaning: after a bond breaks, both molecules cannot bind for this long, so that diffusion can separate
      them (SRSim's answer to geminate recombination); it should satisfy refractory_time x k_off << 1
  - name: k_bond
    type: float
    default: 5.0
    min: 0.0
    range: the published input scripts use fBond = 5.0
    role: thermodynamic
    meaning: spring constant K_d of the harmonic bond potential E_d = K_d (d - d_ij - d_kl)^2
  - name: k_angle
    type: float
    default: 5.0
    min: 0.0
    range: 'published fAngle: 5.0, 15.0, and 0.0 for the scaffold run of figure 6 (angular forces switched
      off)'
    role: thermodynamic
    meaning: spring constant K_a of the harmonic angle potential E_a = K_a (alpha - alpha_ijk)^2
  - name: k_repulsion
    type: float
    default: 5.0
    min: 0.0
    range: the published input scripts use fRepulsion = 5.0
    role: thermodynamic
    meaning: amplitude A of the soft-sphere repulsion E_s = A [1 + cos(pi r / r_c)], r_c = r_i + r_j
  - name: integrator
    type: enum
    default: langevin
    choices:
    - langevin
    - brownian
    role: stochastic
    meaning: 'langevin: the paper''s equation F = F_S - gamma_0 v + sqrt(2 k_B T gamma_0) xi integrated
      with velocities (LAMMPS fix langevin + nve); brownian: its overdamped limit, whose mean square displacement
      is exactly 6 D t'
  - name: orientation
    type: enum
    default: bound-only
    choices:
    - bound-only
    - sampled
    role: stochastic
    meaning: 'how an unbound molecule, which SRSim gives no rotational orientation, is treated: bound-only
      accepts it (SRSim''s own test, which only checks the angles of bonds it already has); sampled accepts
      it with the probability (1 - cos angle_tolerance)/2 that a random orientation offers the component,
      which makes the reactive volume of the paper''s supplement exact'
  - name: well_mixed
    type: bool
    default: false
    role: structural
    meaning: 'run the non-spatial control instead: molecules have no positions and a pair is geometrically
      compatible with probability V_react / V_reactor, which is the well-stirred limit the paper compares
      against BioNetGen in figure 6'
  - name: molecule_types
    type: list
    default: []
    range: 'e.g. [{''name'': ''M'', ''radius'': 0.5, ''mass'': 1.0, ''diffusion'': 1.0, ''sites'': [{''name'':
      ''a'', ''dist'': 1.0, ''theta'': 90, ''phi'': 0}, {''name'': ''b'', ''dist'': 1.0, ''theta'': 90,
      ''phi'': 90, ''states'': [''u'', ''p'']}]}]'
    role: structural
    meaning: 'custom only: elementary species, i.e. the JSON form of SRSim''s .geo file (mass, radius,
      diffusion and the polar coordinates of every component)'
  - name: rules
    type: list
    default: []
    range: e.g. ["'grow' M(a) + M(b) <-> M(a!1).M(b!1) @ 0.5, 0.01", "'phos' A(x~u) -> A(x~p) @ 0.1"]
    role: kinetic
    meaning: 'custom only: BNGL-style rules with microscopic rates; a bare site is free, s!1 is the bond
      of this rule, s!+ bound to anything, s!? untested, s~state an internal state; one rule forms one
      bond, breaks one bond or changes one state, over one or two elementary molecules'
  - name: init
    type: dict
    default: {}
    range: 'e.g. {''M'': 50}'
    role: population
    meaning: 'custom only: how many copies of each molecule type are placed in the box at random positions'
  provides:
  - topology
  - stoichiometry
  - catalysts
  - rate-constants
  - space
  - initial-state
  - mass-conservation
  - sequence-structure-function
  fidelity: reconstructed
  sources:
  - 'Banzhaf, W. & Yamamoto, L. (2015). Artificial Chemistries, section 18.3.3: one paragraph on SRSim
    as a spatial rule-based simulator that ''takes into account the detailed spatial conformation of each
    protein ... and positions each molecule in space such that only neighboring molecules may react with
    each other''.'
  - 'Gruenert, G., Ibrahim, B., Lenser, T., Lohel, M., Hinze, T. & Dittrich, P. (2010). Rule-based spatial
    modeling with diffusing, geometrically constrained molecules. BMC Bioinformatics 11:307 (book [351]),
    doi:10.1186/1471-2105-11-307, open access: https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-11-307
    - the geometry model (fig. 2), the reaction criterion and reactive volumes (fig. 4), the Langevin/bond/angle/soft-sphere
    potentials, the refractory time, and the four applications (scaffold proteins figs. 5-6, microtubules
    fig. 7, spheric self-assembly fig. 8, DNA Sierpinski tiles figs. 9-10).'
  - 'Additional file 2 of the same paper, ''Calculation of Kinetic Parameters from Macroscopic Values'':
    k1_mic = k1_mac; k2_mic = 3 k2_mac / (2 pi eta (3 d0^2 t_dist + t_dist^3)(1 - cos t_ang)^2); the geminate-recombination
    analysis phi = 1 - K2/K1 and the refractory-time corrections [A*] = t_ref k_-1 [C] and [C''''] = [C0]/(1
    + t_ref k_-1). https://media.springernature.com/original/springer-static/esm/art%3A10.1186%2F1471-2105-11-307/MediaObjects/12859_2009_3764_MOESM2_ESM.PDF'
  - 'Additional file 1 (SRSimSrc.zip): the released C++ sources. RuleSys/geometry_definition.cpp computes
    the ideal angle between two components from their polar coordinates (unit vectors, acos of the scalar
    product) and reads DistanceDeviation/AngularDeviation; LammpsModules/compute_reapot_atom.cpp tests
    candidate pairs (distance window (d0 +/- deviDist)^2 and the routine testSiteGeo, which compares the
    angle to every already bound component against its ideal value); LammpsModules/fix_srsim.cpp refracts
    both molecules after a break.'
  - 'Additional file 3 (ExamplesSrc.zip): the published input files. 001.scaffold/spass_w.bngl and spass.geo
    give the scaffold model reproduced here (S with four sites at theta = 0, phi = 0/72/144/216, dist
    8, radius 8, mass 5; A with sites r (theta 0, dist 1, states nP/P) and s (theta 100, dist 1), radius
    1; d1_on 1e-4, d1_off 0.1, d2_on 0.05, d2_off 1e-3, d3_on 1e-5, d3_off 0.9; DistanceDeviation 0.5,
    AngularDeviation 180) and spass_w.in its prefactors (preFactBindR 4.7, preFactBreakR 1e-3, preFactModifyR_1
    1e-3, preFactModifyR_2 4.7), fRepulsion/fBond/fAngle = 5/5/0 and the 100000-step run. Its w_scf.srsim.gdat,
    wo_scf.srsim.gdat, w_scf.gdat and wo_scf.gdat are the data of figure 6: in SRSim 107 phosphorylated
    A with the scaffold against 81 without, in BioNetGen 61 with against 79 without.'
  - 'SRSim manual (docs/Manual.pdf of additional file 1): the .geo/.tgeo file formats and the meaning
    of the fix srsim arguments.'
  decisions:
  - 'The book gives SRSim one paragraph, so everything is reconstructed from the 2010 paper, its supplement,
    its sources and its example models. LAMMPS and BioNetGen are not dependencies: the molecular dynamics,
    the BNGL-style rule language and the rule engine are reimplemented in numpy at the scale that fits
    Chemart''s budget.'
  - 'Rule language: a documented BNGL subset. A rule has one or two elementary molecules per side (SRSim
    executes mono- and bimolecular reactions) and performs exactly one action - form a bond, break a bond,
    or change one internal state - which is SRSim''s set of modifying, binding and breaking rules. Patterns
    constrain one elementary molecule each, not larger subgraphs, so context conditions that span three
    or more molecules (used in the paper''s microtubule and Sierpinski models) cannot be written; exchange
    rules, molecule creation and deletion are not supported by SRSim either.'
  - 'Reaction criterion taken from the released code rather than from the prose: a candidate pair must
    lie in the distance window (d0 - t_dist, d0 + t_dist) with d0 = d_ij + d_kl, and for each partner
    every *already bound* component must keep its ideal angle to the new direction within t_ang (testSiteGeo).
    An unbound molecule has no rotational orientation and passes the angle test; the paper''s supplement
    instead assumes both orientations are fixed when it derives the reactive volume, so orientation=''sampled''
    reproduces that derivation exactly by accepting an unoriented partner with probability (1 - cos t_ang)/2.
    At t_ang = 180 degrees, the value the scaffold model uses, the two agree.'
  - 'Kinetics: a compatible pair reacts with probability 1 - exp(-k dt) (paper, Reaction Model and Kinetics).
    Monomolecular rules are applied per matching molecule with the same expression instead of SRSim''s
    fragmented Gillespie loop over the reactant-pattern vector; within one time step the two are the same
    distribution. Pattern matches are computed once per step and every molecule reacts at most once per
    step, where SRSim recomputes the pattern indices after each reaction.'
  - 'Rate constants in the network: every observed reaction carries the macroscopic constant of the rule
    that fired it, k_mic * V_react for bimolecular rules (a volume per time, so that k [A][B] is the well-stirred
    rate) and k_mic for the others, following the conversion of additional file 2. Reactions produced
    by several rules sum those constants and extras.reaction_rules lists the rules. No symmetry factor
    is applied: as in KaSim''s default convention, each ordered assignment of the rule''s patterns to
    the pair is a separate instance.'
  - 'Space: a periodic box with the minimum-image convention (the scaffold example uses p p p; the other
    examples use reflecting walls, which are not offered). Bonded pairs are excluded from the soft-sphere
    repulsion, as LAMMPS does with special_bonds, and rigid bodies (used for the monomers of the spheric
    self-assembly) are not implemented.'
  - 'The well-mixed control is Chemart''s, not SRSim''s: BioNetGen is not a dependency, so well_mixed=true
    keeps the same rule engine but drops the positions and makes any pair compatible with probability
    V_react / V_reactor, which is exactly the well-stirred limit that additional file 2 derives. This
    is what the figure-6 comparison uses.'
  - 'Scaffold model: the published geometry and rates, with the rate prefactors of the published LAMMPS
    script applied (bind x 4.7, break x 1e-3, monomolecular modify x 1e-3, bimolecular modify x 4.7).
    The published run is 100000 steps; rate_scale compresses it (rate_scale x simulated time is the published
    simulated time), because scaling every rate rescales time. Molecule numbers are scaled down at constant
    concentration. The absolute counts of figure 6 are therefore not reproduced; the effect is: with the
    scaffold the number of phosphorylated A roughly doubles in the spatial simulation and does not increase
    in the well-mixed one.'
  - 'Species identity: the canonical text of the complex (colour refinement, then the least breadth-first
    text over the refined roots, sites grouped by name and bonds numbered by appearance), so isomorphic
    complexes are one species and components with the same name are interchangeable, as in BNGL. The refinement
    is exact for the complexes these models build; texts longer than 120 characters are replaced by formula#hash
    to keep species ids readable.'
  - v1 parameters agent_rules (callable) and geometry (matrix) became the model choice plus the custom
    molecule_types/rules/init, in JSON and in the rule mini-language. The reactor, the tolerances and
    the force constants are exposed with the values the published input files use.
  - 'The network has status ''observed'': it is the set of complex-level reactions that fired, with counts.
    Nothing is enumerated in advance - that is the point of rule-based modelling - so the network is a
    sample of an unbounded species set, and the rule set itself is kept in extras.rules and extras.bngl.'
  phenomena:
  - 'geometry decides the reachable species: paper fig. 3 shows rods from a 180-degree monomer and closed
    squares from a 90-degree one; Chemart''s tests add that with 60 degrees triangles close and with 120
    degrees no small ring can close at all'
  - 'scaffold proteins co-localise their ligands: binding A to the scaffold roughly doubles the phosphorylated
    A in the spatial simulation while the well-mixed simulation of the same rules shows no increase (paper
    figs. 5-6)'
  - in the fast-diffusion limit the spatial simulation reproduces well-stirred mass-action kinetics, with
    the equilibrium constant k2_mic V_react / k_off of the paper's supplement
  - a refractory time suppresses geminate recombination; the molecules waiting it out are t_ref k_off
    [C] per species, as the supplement predicts
  - self-assembly of filaments, rings and larger complexes from a handful of rules, with a combinatorially
    large (unbounded) species set that is never enumerated
  notes: 'The book''s answer to combinatorial explosion: never enumerate. The rule set, the molecule geometries
    and the reaction criterion are kept unflattened in extras (rules, bngl, molecule_types, space.reaction_criterion)
    so that an exporter can hand them to BioNetGen or SRSim itself, while the network reports only what
    actually happened in the reactor.

    '