Table 3 
• Initialize all species and rate constants 
• 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 semiimplicit numerical method 
* Update any rates that have been changed by the continuous reactions and the single discrete reaction 
* Break when userdefined total simulation time is reached 
• end loop 

Adalsteinsson et al. BMC Bioinformatics 2004 5:24 doi:10.1186/14712105524 