Mathematical modelling is an integral part of any study in Systems Biologyand the importance of stochastic effects on the cellular level has been shownfor several systems. Biochemical reaction networks naturally fit thisas the source of stochasticity can be interpreted as the random timing of theunderlying chemical reactions. The dynamics of such a network can be described bya probability distribution that evolves according to the Chemical Master Equation.
In most cases the CME cannot be solved directly butrather has to be estimated by Monte Carlo methods such as the StochasticSimulation Algorithm (SSA). For systems with high timescale separation,either due to variations in magnitude of reaction rateparameters or variation in copy numbers of molecular species, such simulationscan take an impractical amount of computational time.As a remedy one can compute an approximate probability distributionby exploiting such timescale separations, either by regarding reactions asoccuring continuously or using thequasi-stationary assumption to approximate dynamics of fast subnetworks.Existing methods require a manual partitioning of the reaction dynamicsinto discrete or continuous dynamics and a manual identification offast subnetworks. More importantly, if the magnitude of copy numberschange during simulation the identified partitioning and fast subnetworks mightnot be valid anymore, but due to the manual intervention they are fixedin the existing methods.

In this work we propose a computational method to simulate biochemical reaction networkswith multiple timescales.To this end we make use of a concise mathematical framework recently introducedby Kang et al. to simplify the dynamics atdifferent timescales in biochemical networks and we decide the partitioning of species and reactionsusing Linear Programming.This Hybrid scheme is integrated into a simulation algorithm that also allowsthe application of the quasi-stationary approximation for species involved in fast reactionsif their detailed dynamics are not desired by the user.The simulation algorithm is embedded into an adaptive scheme to cope with the problem ofchanging orders of magnitude of copy numbers and shows significant speedups for a number of examplesshowing timescale-separation.