Stochastic simulation: kinetic Monte Carlo and Alchemist

1. From model checking to Monte Carlo

Chapter 11 ended with a stairway: simulation, approximate model checking, exact model checking. This chapter climbs down to the first rung and makes it industrial. Deck 10, "Stochastic Simulation in Alchemist" (Pianini & Viroli), starts where model checking stops:

"Pros [of model checking]: complete exploration of the system, exact verification of property values. Cons: in general extremely costly in terms of memory and time; complexity quickly grows with states; normally only feasible with simple, small systems."

When the state space is too large to explore exhaustively, the Monte Carlo method takes over: "when it's impossible to explore the whole system, find a procedure that randomly explores a part of it, apply it repeatedly, aggregate the result." The name is a trivia with a purpose: "the name is after the famous Casino of Monte Carlo, and refers to the exploration of the probabilities that gamblers can perform by repeatedly playing and recording results."

The crucial precision of the lecture: the procedure can possibly — not compulsorily — be a simulation. The canonical example is geometric: given a rectangle of known area AR containing some irregular figures, sample N points uniformly inside the rectangle, count how many fall inside the figures (say n), and estimate the combined area as AF ≈ (n/N)·AR. "This is not a simulation": no process evolves over time — the same statistical machinery, however, will later drive the simulation of chemical systems, where the "procedure" is a run of the system. The widget below lets you watch the estimate converge.

2. Simulation, reproducibility and time models

With the method in place, the lecture fixes the vocabulary. Simulation is, by the standard definition of Banks et al. (2010), "imitation of the operation of a real-world process or system over time" — and it is not necessarily run on computers: "putting a Formula 1 model into a wind tunnel is a sort of simulation". The model is the imitation itself: "a simplified version of the reality", where simplification is often a requirement because the original process "requires too much time, is not replicable in controlled environments, is too dangerous to replicate, is beyond our technical capacity". One constraint keeps the simplification honest: "elements relevant to the experiment must be retained in the model".

Reproducibility "measures the degree of agreement between two repetitions of the same experiment", and the deck splits it into three notions:

• repeatability — the experiment can be executed multiple times in the same conditions;
• replicability — the experiment can be executed again in slightly different conditions;
• "strict" reproducibility — the computer-science reading: "the ability to obtain the exact same output on the exact same input data".

"A good simulation should be repeatable and replicable. A good computer simulation should be strictly reproducible": randomness must be kept under control — and "parallelism is a source of randomness too!" — a warning that will matter when Alchemist runs large simulations on many cores.

Modelling time

How does a simulation advance time? Two archetypal answers:

• Time-driven — "time is simulated through discrete time slots (ticks); at every tick, the model is updated to reflect the new state; all the changes occurring during the same tick are considered to be simultaneous".
• Event-driven (Discrete Events Simulation, DES) — "events are simulated one by one; for every simulated event, the time is shifted forward; events are strictly ordered: in case two events are scheduled for the same time, one of the two is executed first (and its outcome may influence the remainder of the simulation)".

The distinction is not cosmetic: chemical reactions, message deliveries and node failures do not happen on a shared clock, and the kinetic Monte Carlo of the next sections is fundamentally event-driven — Alchemist is, by design, "a discrete-event simulator: events are forced to be ordered, even if they happen at the same time".

3. Kinetic Monte Carlo: the chemistry problem

The lecture now builds a simulator from scratch: "we have a container with a precise number of molecules that may react with each other. We want to forecast the evolution of the system in the future." Classic chemistry "relax[s] to continuous": differential-equation methods "suppose the concentration of each reactant to be ∈ ℝ", which "is an approximation: you cannot have a quarter of a molecule!", and "these methods are accurate only for a high number of molecules".

What if the system contains "a few thousand molecules"? The Monte Carlo way: "let's start with the system in initial state, let it run and see how it behaves. Repeat." In a real setup this is "very hard to do", so "here it comes the simulation". The system is composed of molecules and reactions of the form

      k
A + B −→ C

where A and B are reactants, k is the reaction-rate constant and C the product. The solution was "first proposed in [Gillespie, 1977] (Gillespie algorithm or Kinetic Monte Carlo)":

1. Compute the propensity ar of each reaction. For a simple bimolecular reaction between distinct species, ar = k[A][B].
2. Execute it, changing the concentrations.
3. Update the propensities which may have changed.

The propensity is the bridge between the model and the math: it is the rate at which the reaction is currently ready to fire — proportional to how often the reactants meet, which is why it is proportional to the product of the concentrations. In the course's own terms, this is the same rate that labels the transitions of a CTMC (chapter 10, section 5): a system of reacting molecules is a CTMC whose state is the vector of molecule counts.

4. The mathematics of the next reaction

To simulate, the algorithm must answer two questions at each step: which reaction fires next, and when. The lecture derives both from a single joint density. For τ ≥ 0, p(τ, μ)dτ is the probability that "no reaction fires in the interval [0, τ) and reaction μ fires in the interval [τ, τ + dτ)". The derivation:

• each reaction j has an exponential time with rate aj, so it "survives up to τ" with probability e−ajτ;
• independence gives P(no reaction fires before τ) = Πj e−ajτ = e−a0τ where a0 = Σj aj;
• μ fires in the next infinitesimal interval with probability aμdτ + o(dτ);
• dividing by dτ and taking the limit: p(τ, μ) = aμ e−a0τ.

From the joint density everything else follows by integration. The probability that μ is the next reaction to fire:

P(next = μ) = ∫₀^∞ a_μ e^(−a_0 τ) dτ = a_μ / a_0

— the propensity, normalised by the total propensity: a weighted coin flip. The waiting time until some reaction fires sums over all reactions, giving the density p(τ) = a0 e−a0τ; its cumulative distribution is F(t) = 1 − e−a0t. Inverting it — "inverse transform sampling" — turns a uniform random number ρ in (0, 1] into an exponential draw:

1 − e^(−a_0 t) = ρ  ⇒  t = −ln(1 − ρ) / a_0 ≡ −ln(ρ) / a_0

The lecture's comment fixes the intuition: "the exponential term comes from the waiting time of the minimum of independent exponential clocks" — each reaction is a Poisson process, and the next firing is the earliest of all the clocks. This is exactly the CTMC view of chapter 10, section 4: the memoryless property is what makes the whole construction valid.

5. The Gillespie algorithm

Everything is in place for the base algorithm. With U(0, 1) denoting the uniform distribution on (0, 1):

1. Set the simulation time T = 0.
2. For each reaction r ∈ R, compute ar, then compute a0 = Σj∈R aj.
3. Draw ρ1 ∼ U(0,1) and select the next reaction μ so that Σj=1μ aj > ρ1 a0 — the cumulative-sum scan.
4. Execute the reaction, changing the concentrations.
5. Draw ρ2 ∼ U(0,1) and advance the simulation time to T = Tprev + (−ln ρ2)/a0.
6. Repeat from step 2.

The data structure of the naive implementation: "choose the next reaction with a cumulative-sum scan over the propensities: sample ρ2a0 and select the first reaction whose cumulative propensity exceeds it (linear time)". Linearity in the number of reactions is the cost that the next two sections attack.

6. Speeding up: dependency graphs and the next-reaction method

Recomputing all propensities at every step is wasteful: "not every reaction affects the speed of every other: for instance, if A + B −k1→ C executes, the propensity of D + E −k2→ F will not be affected." The first optimisation is the dependency graph: "we can improve consistently the performance of the algorithm by keeping in memory which reactions influence which other, and updating only those required" — "a map that connects each reaction to a set of reactions that must be updated".

The deck's example network — five reactions with tangled dependencies:

R1: A + B → C      R2: B + C → D      R3: E + G → A
R4: D + E → E + F   R5: F → D + G

The second optimisation, the next-reaction method (Gibson & Bruck, 2000): "instead of choosing the next reaction probabilistically by propensity, generate a putative time for each reaction; sort the reactions by putative time, and take the first; at each step, for each reaction whose putative time has changed, re-sort the element." The dependency graph is reused to know whose putative time changed. Data structure: "we only need that the first element is the next to be executed — the best solution is a binary heap, which can be accessed in O(1) and sorted in log(n), but with a much smaller average complexity" (in the original work, an "Indexed priority queue").

The third optimisation, random reuse, addresses the cost of generating random numbers — "in a purely chemical simulator, it is often the heaviest task [Gibson and Bruck, 2000]; reducing the number of generated random numbers is key". The trick: "if a dependent reaction did not fire but its propensity changed from ap to ac, update its putative time without drawing a new random number":

τ_c = T + (a_p / a_c) · (τ_p − T)

where T is the current simulation time and τp, τc are the old and new putative times. "This follows from memorylessness: conditioning on survival up to T leaves an exponential residual waiting time, rescaled by the new propensity." The annotation matters: the random reuse is NOT allowed for non-exponential events — a restriction that will matter in section 9, when Alchemist admits events that are not Markovian.

7. Slepoy’s constant-time algorithm

The next-reaction method still costs log n per re-sort. Slepoy's algorithm (Slepoy, Thompson & Plimpton, 2008) pushes selection towards constant time. The idea: "group reactions by propensity magnitude: [pmin, 2pmin), [2pmin, 4pmin), ...; select a group proportionally to its total propensity; then select a reaction inside that group with rejection sampling":

1. Pick a candidate reaction r uniformly from the group.
2. Draw u ∼ U(0,1); accept r if u < ar/amax, where amax bounds all propensities in the group.
3. If rejected, draw a new candidate and try again.

The cost analysis: "on average, rejection needs ≤ 2 trials (propensities in a group differ by at most a factor 2); if the number of groups G stays bounded, reaction selection is O(1)". The deck is careful about the assumptions: "the O(1) claim assumes the propensity range stays bounded; constant-time updates also require bounded coupling: each fired reaction should affect only O(1) other propensities." The lecture's lesson for the engineer: constant-time selection is achievable, but it rests on structure — bounded propensity ratios inside groups, bounded fan-out in the dependency graph — which is precisely the structure Alchemist's spatial dependency graph (section 9) tries to guarantee.

8. Alchemist: from chemistry to pervasive computing

The lecture now turns to the tool the course uses for large-scale simulation, Alchemist. The background: "pervasive computing scenarios are normally simulated by means of 'agent-based simulators' (ABS) [Wooldridge and Jennings, 1995]. ABS are extremely flexible, but they lack performance: it's the price to pay for being able to simulate a very wide spectrum of situations." The key observation that reframes the problem:

"Many pervasive computing scenarios can be modelled as mobile multi-compartmented chemical systems, where molecules are pieces of data (equivalent to a network of Petri Nets)."

This is the link with chapter 9 — a network of communicating computational devices is structurally a Petri Net with places as compartments, tokens as data and transitions as reactions. "A whole literature exists on how to make very fast kinetic Monte Carlo algorithms", so instead of optimising an ABS at the simulation level, the research question is: "can we take a kinetic Monte Carlo and extend it until it supports all the abstractions we need?"

The target scenarios: self-organising systems, pervasive computing systems, crowds of people, large-scale situated systems, smart mobility, crowd detection and steering, sensor networks, computational biology, aggregate programming. The requirements that a simulator must then meet:

• multiple compartments (from now on: nodes);
• molecules can be different data types;
• node mobility;
• non-Markovian events;
• a more flexible concept of reaction;
• high performance.

These requirements rewrite the abstract model. The deck's vocabulary is precise and maps one-to-one onto the KMC ingredients: the Environment is "a Riemannian manifold where nodes live"; a Node is "a container of reactions and molecules situated in the environment"; a Molecule is "a token representing a chunk of data (think of it as a pointer)"; the Concentration is "the actual data associated with a molecule"; a Reaction is "a proactive behaviour"; and a Linking Rule is "a function of the environment that decides whether or not two nodes are connected".

A reaction, in turn, is tripartite: conditions ("node contains something", "number of neighbors < 3", any other condition about the environment), a probability distribution ("rate equation: how conditions influence the execution speed"), and actions ("change concentration of something", "move a node towards...", any other action on the environment). The propensity of section 3 is the rate equation in disguise: the conditions compute the concentration factors, the distribution turns them into a firing rate.

9. The engine: spatial dependencies and non-Markovian events

Moving from one container to many changes the dependency graph. With multiple compartments: "there are more reactions: each node has its 'copy'; a reaction may affect the propensities locally, in the neighborhood, or globally; the fewer are the bindings between reactions, the higher is the efficiency of a dependency graph". The lecture's solution is the spatial dependency graph, built on three contextual levels — local, neighborhood, global:

• each reaction is assigned an input context — which parts of the environment it reads to compute its propensity;
• each reaction is assigned an output context — which parts of the environment it modifies;
• a reaction r1 may influence r2 if one of the following holds:
  • r1 and r2 belong to the same compartment;
  • r1's output context is global;
  • r2's input context is global;
  • r1's output context is neighborhood and r2 belongs to the neighborhood;
  • r2's input context is neighborhood and r1 belongs to the neighborhood;
  • both r1's output and r2's input contexts are neighborhood, and there is a compartment shared by the two neighborhoods.

The rule set is a precise engineering answer to the three challenges of multi-compartment simulation: "who decides if two compartments are communicating?" (the linking rule), "how to model a molecule moving towards a new node?" (a reaction whose action moves data — or a node), "how does the dependency graph change?" (only along the contexts above).

Non-Markovian events

The second engine extension is about time. Example: "every second, an external device injects some quantity of molecules within a compartment — this event happens precisely every second: it is not a Poisson process! Its probability distribution is a δ-Dirac comb." The basic Gillespie algorithm "is hard to modify to support such events. The main reason is that the choice is not made depending on time, but on propensity, which is an entity strictly bound to the Markovian model." The next-reaction algorithm, instead, "uses putative times: this makes it able to simulate events independently from their distribution, since we just need to correctly estimate the next time of occurrence." With the restriction already met in section 6: "the random reuse is NOT allowed for non-exponential events."

10. The SAPERE incarnation and YAML simulations

Alchemist was "initially developed within the SAPERE EU Project; at the model level, it captures the required abstractions of a SAPERE system". In the SAPERE incarnation, "the concentration is defined as 'list of tuples matching a tuple template'. Basically, in this configuration Alchemist does not simulate a simple collection of intercommunicating compartments, but a network of (possibly mobile) programmable tuple spaces." This is "one of four incarnations, along with two aggregate programming incarnations (Protelis and Scafi) and a sketched biochemical implementation" — the last two are the subject of chapter 14.

Simulations are written in YAML: "JSON superset, compact, human-readable, supports anchoring (referencing)"; "the same syntax can be used for any incarnation"; there is "support for running batches, no intermediate compilation, no large files involved". The mini-tutorial of the deck builds a simulation from the ground up.

incarnation sapere

network-model
  type ConnectWithinDistance   # a LinkingRule implementation
  parameters [5]               # connect nodes within distance 5

deployments
  - type Point                 # a Displacement implementation
    parameters [0, 0]
  - type Point
    parameters [0, 1]

incarnation sapere

network-model
  type ConnectWithinDistance
  parameters [0.5]

deployments
  type Grid
  parameters [-5, -5, 5, 5, 0.25, 0.25, 0.1, 0.1]
  contents
    - molecule hello              # everywhere
    - in                          # restrict the area...
        type Rectangle
        parameters [-1, -1, 2, 2] # ...to this rectangle
      molecule token

programs
  -
    - time-distribution 1         # Markovian rate 1
      program >                   # '>' begins a multiline string
        { token } --> { firing }
    - program "{ firing } --> +{ token }"   # ASAP reaction

incarnation sapere
network-model
  type ConnectWithinDistance
  parameters [0.5]

_send & send                       # anchor: a named program list
  - time-distribution 1
    program >
      { token } --> { token } * { token }    # send a copy to a neighbor
  - program >
      { token }{ token } --> { token }       # merge duplicates

deployments
  type Grid
  parameters [-5, -5, 5, 5, 0.25, 0.25, 0.1, 0.1]
  contents
    in
      type Rectangle
      parameters [-0.5, -0.5, 1, 1]
    molecule token
  programs *send

_grad & grad
  - time-distribution 0.1
    program "{ source } --> { source } { gradient, 0 }"
  - time-distribution 1
    program "{ gradient, N } --> { gradient, N } *{ gradient, N +# D }"
  - program >
      { gradient, N }{ gradient, def N2 >= N } --> { gradient, N }
  - time-distribution 0.1
    program "{ gradient, N } --> { gradient, N + 1 }"
  - program "{ gradient, def N > 30 } -->"

programs
  - time-distribution
      type DiracComb
      parameters [0.5]
    program "{ token, N, L } --> { token, N, L } *{ token, N +# D, L add [#NODE;] }"
  - program "{ token, N, L }{ token, def N2 >= N, L2 } --> { token, N, L }"

11. Lab: operational steps and R&D directions

Deck 10 has no separate lab deck: the tutorial is the lab, and the lecture's UI notes close the operational picture. The Alchemist GUI keyboard shortcuts: P pause/play, L toggles link painting, R realtime mode ("tries to sync the simulation with the real time, always ensuring at least 25fps"), →/← speed up / slow down (fewer / more UI updates), M marker for the node closest to the mouse, S select mode, O manual move of selected nodes.

Operational steps

• Step 1 — get Alchemist running. Download the latest Alchemist distribution (see alchemistsimulator.github.io); write the minimal incarnation sapere simulation from section 10 and open it in the GUI. Add deployments (Point, Circle, Grid), then contents and programs, and watch the reaction events fire in the trace panel. Verify the UI shortcuts above.
• Step 2 — reproduce the tutorial. Build the diffusion program and the grad program exactly as in the deck; seed source molecules in a rectangle and observe the gradient field spreading. Then add the DiracComb injection and check that injections happen exactly every 0.5 time units — the non-Markovian event of section 9.
• Step 3 — parametrise. Move the simulation to the variables section (rate, size, sourceStart, sourceSize) and run a batch sweep; compare how the spreading time changes with size. This is the operational counterpart of the PRISM experiments of chapter 11, section 6 — parametric analysis, by simulation instead of model checking.

R&D tasks

Test your knowledge