Table 3

• Compute all reaction rates

• Loop:

   * Set μ = sum of rates for the discrete reactions

   * if (p _t = μΔt > ε), use Gillespie algorithm:

* R = a uniform random number in (0,1)

* Set timeStep = -log(R)/μ

* Find which reaction occurred, update the species involved

   * else, use small Δt approximation:

* R = a uniform random number in [0,1]

* timeStep = continuousTimeStep

* if (R <p _t = μ × timeStep), discrete transition has occurred:

• Determine which discrete transition occurred:

• Find the first value of k for which

• If , the forward reaction occurred, otherwise the backward reaction occurred

* else, no discrete transition:

• No discrete reaction occurs, update is entirely due to continuous reactions (below)

* end if (small Δt method, determination if discrete transition occurred)

   * end if (selection of Gillespie or small Δt method for discrete reactions)

   * Update the continuous species using the Langevin equation, with step size timeStep (where timeStep is either equal to continuousTimeStep or to the step size found by the Gillespie algorithm), using a semi-implicit numerical method

   * Update any rates that have been changed by the continuous reactions and the single discrete reaction

   * Break when user-defined total simulation time is reached

• end loop