Given the number of the different species X(t 0 ) present at a time t 0 the next reaction to occur and the time ∆t 0 before it occurs is computed as follows. Each possible reaction j has an associated probability to occur equal to p j (t 0 ) = a j (X(t 0 ))/a(X(t 0 )) where a j (X(t 0 )) is the propensity function of the reaction j given the state of the system X(t 0 ) and a(X(t 0 )) = i a i (X(t 0 )) is the sum of all propensity functions at time t 0 . The time ∆t 0 at which the next reaction will occur is randomly sampled from an exponential distribution p(∆t 0 |X(t 0 )) = a(X(t 0 )) exp[-a(X(t 0 ))∆t 0 ]. In the "atomic-features" framework, the atoms involved in the randomly-selected reaction must also be picked. Each species with the fingerprint associated with the selected reaction has equal weight to undergo the reaction. The bonds of the system are modified according to the selected reaction and atoms, the number of the different species X(t 0 ) is updated to X(t 0 +∆t 0 ). This entire process is then repeated until a simulation time t is reached.