HSimulator : Hybrid Stochastic / Deterministic Simulation of Biochemical Reaction Networks

HSimulator is a multithread simulator for mass-action biochemical reaction systems placed in a well-mixed environment. HSimulator provides optimized implementation of a set of widespread state-of-the-art stochastic, deterministic, and hybrid simulation strategies including the first publicly available implementation of the Hybrid Rejection-based Stochastic Simulation Algorithm (HRSSA). HRSSA, the fastest hybrid algorithm to date, allows for an efficient simulation of the models while ensuring the exact simulation of a subset of the reaction network modeling slow reactions. Benchmarks show that HSimulator is often considerably faster than the other considered simulators. The software, running on Java v6.0 or higher, offers a simulation GUI for modeling and visually exploring biological processes and a Javadoc-documented Java library to support the development of custom applications. HSimulator is released under the COSBI Shared Source license agreement (COSBI-SSLA).

We introduce HSimulator, a Java hybrid stochastic/ deterministic simulator for mass-action biochemical reaction systems placed in a well-mixed environment, where position and speed of molecular species are randomized and therefore they do not affect reaction executions.Species  1 , . . .,   are represented in terms of abundances (number of molecules) and reactions  1 , . . .,   are defined as where the species on the left of the arrow are reactants and the ones on the right are products.The stoichiometric coefficients V  and V   indicate how many molecules of reactant are consumed and how many molecules of product are produced, respectively.The constant   at the end of the reaction is the stochastic reaction constant introduced by Gillespie [8] for computing reaction firing according to the definition of massaction propensities.
From the milestone work of Gillespie [8], where the Direct Method (DM) has been defined, exact stochastic simulation is the most accurate simulation approach.Its drawback is the high computational complexity arising from the need of separately simulating each reaction firing.Several algorithms have been introduced to decrease simulation runtime, such as the Next Reaction Method (NRM, [9]) and later the Rejectionbased Stochastic Simulation Algorithm (RSSA, [10]).The latter is tailored for complex biochemical reactions with timeconsuming propensity functions and it constitutes the state of the art of exact stochastic simulation.
Despite the improvements, the problem of simulation complexity remains.In fact, the investigation of complex diseases or disorders is often concerned with extended biomolecular networks involving genes, proteins, metabolites, and signal transduction cascades.This justifies the introduction of approximate simulation strategies, which sacrifice accuracy to decrease runtime.Such strategies progressively reduce the stochasticity of the dynamics [11,12] until reaching a deterministic simulation, where the model is translated to a 2 Complexity system of ordinary differential equations (ODEs, [7,13]) and the dynamics provides an averaged behavior.
A drawback of an approximate simulation is that all reactions in the model are simulated according to the same approximate strategy.This is an issue when exact simulation is required at least for a small part of the reaction network, for example, when slow reactions are considered to model rare stochastic events.In such a case, a hybrid simulation strategy is able to partition reactions into subsets that are simulated by different algorithms [7,14,15].Hybrid simulation therefore needs to synchronize the progress of different simulation strategies in order to guarantee the exactness of the simulation of slow reactions.This requires considering time-varying transition propensities, which in turn require the computation of integrals to calculate the exact firing time of slow reactions [7,14,[16][17][18][19][20] (see also Section 2.5 and (19)).This step is computationally demanding and practical implementation necessarily approximates the integral computation by a numerical method that affects the accuracy of the simulation of slow reactions.
The Hybrid Rejection-based Stochastic Simulation Algorithm (HRSSA) has been recently introduced in [21] to avoid approximations in simulating slow reactions.HRSSA exploits some computational advantages introduced in RSSA to exactly simulate slow reactions and to apply an efficient and accurate dynamic partitioning of reactions, which constitute the two most significant bottlenecks of hybrid simulation.
HSimulator herein presented fills a gap in the current literature of stochastic and hybrid simulation by providing a suite of published state-of-the-art simulation algorithms including, in the same package, the exact algorithm RSSA and the first publicly available implementation of its hybrid version HRSSA.The benchmarks in the paper show that the HSimulator implementation is often faster than the state-ofthe-art simulator COPASI [22].HSimulator is released under the COSBI Shared Source license agreement (COSBI-SSLA) and it is available for download at the COSBI website at http://www.cosbi.eu/research/prototypes/hsimulator.

Materials and Methods
HSimulator provides an implementation of five state-of-theart simulation strategies covering exact stochastic simulation (DM and RSSA), deterministic simulation (Forward Euler and the Runge-Kutta-Fehlberg RK45 adaptive algorithm), and hybrid simulation (HRSSA).For deterministic and hybrid simulations, the simulator automatically translates the mass-action reaction network into a set of ordinary differential equations (ODEs, see Section 2.3).In the following some insights about the implemented simulation algorithms are provided as well as their pseudocodes, which can be used as a reference to better understand the functionalities implemented in the simulator (for additional details, the reader is invited to refer to the user guide of HSimulator available at http://www.cosbi.eu/research/prototypes/hsimulator).Readers already familiar with the topic can skip this section and going directly to Section 3.

Direct Method (DM).
The Direct Method constitutes the first practical implementation of an exact stochastic simulator [8].The simulation of an exact algorithm is based on the concept of reaction probability density function (pdf): which provides the probability (,  | x, ) that the next reaction to apply will be   in the infinitesimal time interval [ + ;  +  + ], given the system state x at time .
The state of the system at time  is represented by the column vector where   () is the number of molecules of species   in the system at time .In the following we will often write x for X() and   for   () to improve the readability of formulas.
The term   (x) in ( 2) is the propensity of   in the state x, while the total propensity  0 (x) is the sum of the propensities of all reactions: The reaction propensities   (x) are computed from the stochastic reaction rate   : where ℎ  (x) is the number of distinct reactant combinations for   in the state x [8].For standard mass-action kinetics, as considered here, it is where V  are the stoichiometric coefficients of   according to (1).Equation (2) shows that the next reaction   fires with a discrete probability   (x)/ 0 (x) and its firing time  has an exponential distribution with rate  0 (x).Gillespie introduced the Direct Method (DM) for exactly sampling the pdf (,  | x, ) by applying the inverse transformation: where  1 and  2 are random numbers drawn from a uniform distribution (0, 1).The pseudocode of the DM algorithm is in Algorithm 1.
Input: a reaction system with  species and  reactions as the one in (

The Rejection-Based Stochastic Simulation Algorithm (RSSA).
The Rejection-based Stochastic Simulation Algorithm (RSSA) is a novel exact stochastic simulation algorithm introduced in [10] and further improved in [7,[24][25][26].RSSA constitutes the state of the art of exact stochastic simulation tailored for complex biochemical reactions with time-consuming propensity functions.The pseudocode of RSSA is provided in Algorithm 2.
RSSA computes for each reaction a lower bound and an upper bound of the propensity, which encompass all possible Complexity values of reaction propensities over the fluctuation interval.Such bounds are derived by defining a fluctuation interval X = [x; x] for the system state x (step (2)), where x = ⌊x − x⌋ , and 0 <  < 1 is a simulation parameter whose value is usually between 0.1 and 0.2 (10-20% of the system state).Because reaction propensities with mass-action kinetics are monotonic functions, the propensity bounds of reaction   , for  = 1, . . ., , correspond to the interval [  (x),   (x)] (step (3)).These propensity bounds are recomputed only when the current system state is not anymore inside its fluctuation interval (x ∉ X).
The selection of reaction firing in RSSA is composed of two steps.First, it selects a candidate reaction   with probability   (x)/ 0 (x) (step (7)), where  0 (x) = ∑  =1   (x).The candidate reaction is then validated according to a rejectionbased procedure that implements the toss of a biased coin with success probability   (x)/  (x) (steps (8)-( 14)).To do this, the algorithm generates a random number  in (0, 1) and moves in one of the following three cases: (1) If  ≤   (x)/  (x), then   is accepted without computing its propensity   (x) because   (x)/  (x) ≤   (x)/  (x) (quick accept); After a reaction is accepted to fire, the state is updated and the algorithm moves to the next simulation iteration without updating the propensity bounds.Only for uncommon cases when state x moves out of its fluctuation interval does RSSA have to define a new fluctuation interval and update the propensity bounds.Numerical experiments show that propensity updates in RSSA are infrequent; hence their impact is very low during the computation [10].This advantage is mitigated by the penalty paid to prepare a potential advancement of the system for a candidate reaction that is finally rejected to fire, but the overall performance is still considerably better than DM and other stochastic simulation algorithms.

The Forward Euler Algorithm.
The forward Euler algorithm is the simplest example of deterministic simulation algorithm.Here the mass-action model is translated to a system of ordinary differential equations (ODEs, [7,13]) and the dynamics provides an averaged behavior of the system.
Consider a biochemical reaction system with  1 , . . .,   species,  1 , . . .,   reactions defined according to (1).The mass-action deterministic rate constant   of   can be obtained from the stochastic one   as where   indicates Avogadro's number,  is the volume where the reaction is taking place, and Order  is the overall order of reaction   , which in mass-action is defined as the sum of the stoichiometric coefficients of reaction reactants.Finally, the ODE describing the evolution in terms of molar concentrations of species   ,  = 1, . . ., , is where squared brackets are used to indicate molar concentration of species ([  ] =   /(  ),  = 1, . . ., ).Once the model has been translated to a system of ODEs a first-order approximation of the next system state can be computed by where ℎ is a user-defined discretization stepsize.The pseudocode of the Forward Euler algorithm is provided in Algorithm 3.

Runge-Kutta-Fehlberg Algorithm (RK45).
The forward Euler algorithm introduced in the previous section is a firstorder numerical method provided for didactic purposes.Several numerical methods have been introduced to increase the accuracy of deterministic simulations and to decrease simulation runtime.A popular algorithm is the Runge-Kutta-Fehlberg algorithm (RK45).This algorithm represents the standard choice of several numerical integrators when the model does not exhibit stiffness.It is also often used in hybrid simulation strategies for the simulation of the fast part of the network.The reader can find a comprehensive guide to numerical methods of ODEs in [7,13].
The Runge-Kutta Fehlberg method of fourth-order, also known as the RK45 method, updates the system state by Input: a system of ODEs [X]/ = F(, [X]) corresponding to a biochemical reaction system, the initial state [X( 0 )] of the system with species concentrations at time 0, the last time instant  end to be simulated and the discretization stepsize ℎ.
Given a system of ODEs as in (11) and a discretization stepsize ℎ, the values of  1 , . . .,  6 are shared between ( 13) and ( 14) and they can be computed as The fourth-order version of the algorithm provided in ( 13) is used to compute the dynamics of the system.The fifth-order scheme of ( 14), instead, is used to estimate the maximum local truncation error introduced in the simulated step: where max provides the maximum value along the vector of truncation errors.The error estimate Δ +1 is then compared to the error threshold   specified by the user.When Δ +1 ≤   , the local truncation error is assumed to be smaller than the threshold, the state [X +1 ] is accepted, and the algorithm moves one step forward.Otherwise, the new state is not accepted and the next state is evaluated again using a different (smaller) value of ℎ.In both cases, the value of ℎ is updated as When the estimations [X +1 ] and [ X+1 ] agree to more significant digits than required, then  becomes greater than 1 and ℎ is increased.Equation ( 18) is derived from the general formula  = (  /Δ +1 ) 1/ , which defines how to update the value of ℎ of an adaptive one-step numerical method of order , by considering the error estimate Δ +1 and the user-defined threshold   .The additional multiplicative factor of 0.84 is an empirical number commonly added in RK45 implementation to reduce the variability of ℎ, because very high values of the stepsize increase the probability of repeating the next computed step.The implementation of the algorithm is in Algorithm 4. The value of ℎ is updated at each step starting from an userprovided initial value ℎ 0 .The next computed state of the system is accepted only when the estimate of the local truncation error Δ remains below the user-provided threshold   .Steps (16)-( 21) are additional steps added to the implementation in order to avoid very large modifications of ℎ in a single step.
Even if the computation of a simulation iteration of the RK45 algorithm is more computational demanding than other nonadaptive Runge-Kutta implementation, the possibility of changing the value of ℎ often dramatically decreases the simulation runtime.

Hybrid Rejection-Based Stochastic Simulation Algorithm (HRSSA).
A drawback of approximate simulation is that all the model reactions are simulated according to the same approximate strategy.This is an issue when exact simulation is required at least for a small part of the reaction network, for example, when slow reactions are considered to model rare stochastic events.In such a case, a hybrid simulation strategy can be applied, which divides reactions into subsets that are simulated by different strategies at the same time [7,14,15].
An issue of this approach is the synchronization between the employed simulation strategies in order to guarantee the Input: a system of ODEs [X]/ = F(, [X]) corresponding to a biochemical reaction system, the initial state [X( 0 )] of the system with species concentrations at time 0, the last time instant  end to be simulated, an initial value for the discretization stepsize ℎ 0 and a threshold for the maximum local truncation error   .Output: a time series of states [X()],  ∈ [ 0 ;  end ], providing the dynamics of the system in terms of molar concentrations.

Pseudocode:
(0) initialize time  = 0, state [X] = [X( 0 )] and discretization stepsize ℎ = ℎ 0 ; (1) while  <  end do exactness of the simulation of slow reactions.This is a fundamental aspect of hybrid simulation relying on the definition of the probability density function (pdf) of slow reactions.In fact, even though fast reactions can be safely simulated across slow reaction events, it is not always the case that slow reactions can be simulated regardless of what fast reactions are changing in the system.Actually the reaction propensity   of a slow reaction   ∈ R slow may depend on species whose quantities are changed by fast reactions.For this reason, the reaction probability density function of (2) is not suitable to simulate the reaction subnetwork made of only the slow reactions of the system and it has to be extended to consider time-varying transition propensities [7,14], which account for the fact that reaction propensities of slow reactions may change over time by the simulation of fast reactions.
According to [7,14,[16][17][18][19][20], the pdf of the next firing of a slow reaction   ∈ R slow finally becomes where The firing time  of the next slow reaction   is thus obtained by solving the equation where  is a random number from (0, 1).Solving ( 21) is computationally challenging because the state () is changed by fast reactions during time interval [,  + ].
Therefore, the hybrid simulation has to evaluate the integral simultaneously with the simulation of fast reactions in order to correctly generate the next slow reaction event.Moreover, in practical implementation this step has to be necessarily approximated by a numerical method relying on an error threshold that introduces an additional approximation in the simulation of slow reactions.This step is so critical that the standard choice of several implementation available in literature, such as the hybrid algorithms implemented in COPASI [22], intentionally avoids the computation of the integral during the simulation by accepting some approximation also in the simulation of slow reactions.The hybrid algorithm HRSSA is a novel hybrid simulation algorithm introduced in [21] to solve this problem.HRSSA is built on top of RSSA introduced in Section 2.2.In particular, the RSSA concept of fluctuation interval of the system Input: The same as RSSA (Algorithm 2) together with the time granularity  fast used for the (approximate) simulation of fast reactions; the minimum amount  ∈ N of molecules that has to be available for fast reactions; the minimum number of times  ∈ N that a fast reaction has to be applied, in average, within the time range of size  fast .Output: the sets R slow and R fast of slow and fast reactions, respectively.Pseudocode: (0) R slow fl 0; R fast fl 0; ( state is widely used to provide an important computational advantage.HRSSA synchronizes the simulation of slow and fast reactions by avoiding the computation of the integral of ( 21) and applies an accurate dynamic partitioning of reactions without updating reaction classification at each simulation step.
HRSSA updates reaction partitioning only when the current system state does not fit anymore in its fluctuation interval (see (8)).This permits reducing the computational overhead without losing the accuracy of the classification.HRSSA considers both reaction propensities and the number of transformed molecules for reaction partitioning, by replacing real propensities with their lower bounds.The adoption of bounds instead of real propensities does not affect the accuracy of the classification.Conversely, the usage of lower bounds imposes more stringent constraints that tend to increase the number of reactions that are classified as slow (and therefore simulated without approximations).The pseudocode of the dynamic reaction partitioning of HRSSA is in Algorithm 5.The if clauses of steps (2) and (4) implement the two conditions on reaction propensities and number of transformed molecules, respectively.
After reaction partitioning in the two sets of slow and fast reactions, HRSSA computes the sum of upper propensity bounds  slow 0 (x) of slow reactions as defined in (20).The firing time of a candidate slow reaction is then computed as where  is a random number in (0, 1).Under the hypothesis that the system state will remain inside its fluctuation interval in [,  + ], we can consider  slow 0 (x) not dependent on time over [, +] and this allows us to simulate fast reactions over this interval without taking any side effect on the application of slow reactions into account.This is because (22) remains valid, regardless of the action of fast reactions, as long as the current system state respects its bounds.
After the simulation of fast reactions for the time interval [,  + ], a slow reaction is chosen and validated to fire by a rejection test, according to the RSSA simulation strategy.To guarantee the exact simulation of slow reactions (the proof is in [21]), the simulation of fast reactions is required to happen in a feasible system state.Therefore, every time the system state exits from its bounds the simulation is stopped and the fluctuation interval is updated.The pseudocode of HRSSA is provided in Algorithm 6.

Results and Discussion
HSimulator is a cross-platform simulation software written in Java, offering a suite of state-of-the-art stochastic, deterministic, and hybrid simulation strategies for massaction well-stirred biochemical reaction systems.Multicompartmental modeling is natively supported allowing the definition of nested reaction volumes via a customary text-based representation of reactions and compartments.The multithreading implementation allows scaling very well when computing averaged model dynamics obtained by running multiple stochastic simulations in parallel.The software has been designed to run in three user-scenarios: (i) via command-line, for example, as batch jobs or as part of wider modeling terminal scripts; (ii) programmatically embedded via its API in custom applications [27] or within MATLAB/Octave/R/Mathematica projects; (iii) as a self-standing simulation environment with a graphical user interface (GUI) for biomedical system modeling.For the three usage scenarios extensive documentation is available at the software web page (http://www.cosbi.eu/research/prototypes/hsimulator.).The GUI assists the modeler while defining the reactions with syntax highlighting of the typed entities and quick graphical access to simulation parameters.The simulation environment Input: The same as RSSA (Algorithm 2) together with the time granularity  fast used for the (approximate) simulation of fast reactions; the minimum amount  ∈ N of molecules that has to be available for fast reactions; the minimum number of times  ∈ N that a fast reaction has to be applied, in average, within the time range of size  fast .Output: a time series of states X(),  ∈ [ 0 ;  end ], providing the dynamics of the system.Pseudocode: (0)  fl  0 ; (1) while  <  end do (2) Define the fluctuation interval X = [x; x] according to ( 8 [21] for details).In HSimulator, step (9) is implemented by a deterministic Runge-Kutta numerical method.is completed by an interactive viewport to explore and understand the modeled system dynamics (see Figure 1) which can be then exported as Excel spreadsheets.
To evaluate the performance of the HSimulator implementation, a set of benchmarks have been carried out considering the state-of-the-art simulator COPASI [22], version 4.19 (build 140).All the simulations have been run in the same conditions on a 64 bit macOS Sierra MacBook Pro, with 16 GB of RAM.Six models have been simulated according to 6 different algorithms covering all the simulation strategies (stochastic, deterministic, and hybrid).Namely, we considered the Direct Method (implemented in both COPASI and HSimulator) and RSSA (only implemented in HSimulator) for exact stochastic simulation, LSODA (implemented in COPASI) and RK45 (implemented in HSimulator) for deterministic simulation, and Hybrid Runge-Kutta (implemented in COPASI) and HRSSA (implemented in HSimulator) for hybrid simulation.Benchmarks have been computed by averaging the runtime of 20 simulations.The simulation parameters for each algorithm are in Table 1 and they have been chosen to preserve as much as possible the reliability of comparisons.
The Hybrid Runge-Kutta algorithm provided by COPASI combines two fast simulation strategies available in literature to simulate slow and fast reactions (NRM combined with the 4th order Runge-Kutta method [7,13]).Reaction partitioning is computed dynamically during the simulation according to a threshold which provides the maximum number of reactant molecules of a slow reaction.Moreover, the simulation of slow and fast reactions is synchronized without computing the integral of ( 21) by approximating reaction probabilities of slow reactions as constant during one stochastic step.This solution introduces an error in the simulation of slow reactions, but makes the computation faster.HRSSA resulted to be faster than this algorithm in all the simulation cases we considered, even if HRSSA does not apply such approximation in simulating slow reactions.
We note that the latest version of COPASI (version 4.19 build 140) considered in our benchmarks provides other two implementation cases of hybrid simulation algorithms (Hybrid LSODA and Hybrid RK45).Here we considered Hybrid Runge-Kutta because its implementation is closer to the one of HRSSA.However, we obtained very similar simulation runtimes also by running our benchmarks with Figure 1: The HSimulator graphical user interface allows easily defining, even complex, biological models in terms of customary text-based reactions.The appropriate simulation strategy can be selected and parameterised as well according to the model, although the HRSSA method offers the best hybrid simulator.The figure shows the Oregonator model [8] simulated by the HRSSA algorithm [21].The plot shows variable abundances in logarithmic scale to illustrate the switch between the deterministic simulation of fast reactions (high abundances) and the exact stochastic simulation of slow reactions (low abundances).Hybrid LSODA (data not shown).Hybrid RK45 has been not considered here, because this simulation strategy uses a static reaction partitioning whereas HSimulator and the compared methods use a dynamical approach allowing for a fair benchmark between strategies.We also have to acknowledge that COPASI is not the only simulation software available in the literature.In the present work we considered COPASI due to its high popularity and to allow, as much as possible, a fair comparison with HSimulator.In fact, COPASI provides a way of specifying the biochemical system that is very close to the one of HSimulator.Snoopy [28] is another promising simulation software that provides different simulation strategies for hybrid models, including a reinterpretation of HRSSA.Simulation benchmarks are provided in Table 2. To allow the reproducibility of the presented results, all the considered models are provided with HSimulator both implemented for HSimulator and for COPASI.Simulation benchmarks comprehend two biological models (the MAPK cascade and the Gemcitabine mechanism of action) to test simulation strategies in real modeling applications plus two theoretical models (the fully connected model and the multiscale model) considered with two different parameterisations each, to evaluate the performance of simulation algorithms under specific conditions.The fully connected model has been considered to evaluate the performance of simulation algorithms when a mass-action biochemical reaction network of reactions that are all slow (MCM  = 20) or fast (MCM  = 2000) is simulated.Conversely, the multiscale model allowed testing the intermediate condition, where the biochemical network can be divided into two subnetworks working at different time scales.This scenario has been specifically considered to test hybrid simulation strategies.The biochemical network of the multiscale model has been generated two times to test the scalability of simulation algorithms (MSM (10, 50), network of 60 species and 140 reactions; MSM (20, 100), network of 120 species and 480 reactions, see Section 3.4  [29].The MAPK cascade is at the heart of a molecular-signaling network that governs the growth, proliferation, differentiation, and survival of many cell types.Moreover, the MAPK pathway is deregulated in various diseases, ranging from cancer to immunological, inflammatory, and degenerative syndromes, and thus represents also an important drug target [30].
The MAPK cascade is a series of three protein kinases that are responsible for cell response to growth factors.The signal 1 activates MAPKKK by phosphorylation, which in turn activates MAPKK.Once activated, MAPKK activates MAPK.When 1 is added to the system, the output of activated MAPK increases rapidly.By removing the signal 1, the output level of activated MAPK reverts back to zero.Reverse reactions are triggered by the signal 2 and by the MAPK phosphatases KKpase and Kpase.The reaction-based representation is A simulation benchmark of the MAPK cascade is in Table 2 (first column).In the computed benchmark, the model has been simulated for 150 time units with all stochastic reaction constants equal to 0.001.The initial abundances have been set to 1 = 2000, 2 = 2000, K = 200000, KK = 200000, and KKK = 20000.All the other species are initialized with 0 molecules.

The Gemcitabine Model.
Gemcitabine is a drug for nonsmall-cell lung cancer, pancreatic cancer, bladder cancer, and breast cancer.Here we will consider a simplified model of its mechanism of action according to [23].
Gemcitabine (dFdC, see Box 1) is transported from plasma into the cell through the cell membrane (dFdC out → dFdC).Gemcitabine can be deaminated by CDA in the cytoplasm and in the extracellular environment leading to dFdU.Both dFdC and dFdU can be phosphorylated by dCK.Monophosphorylated Gemcitabine can be deaminated by dCMPD, whereas dCMPD is inhibited by dFdC TP .Alternatively, it is further phosphorylated.The Gemcitabine triphosphates dFdC TP and dFdU TP compete with the natural nucleoside triphosphate dCTP for incorporation into nascent DNA chain and inhibit DNA synthesis, thus blocking cell proliferation in the early DNA synthesis phase.
A simulation benchmark of the Gemcitabine model is in Table 2 (second column, simulation length 12 time units).In the computed benchmark, stochastic reaction constants are [23]: The fully connected model is useful to evaluate the scalability of algorithms (variation in runtime with growing values of ).Moreover, if we fix initial values of the species and stochastic reaction constants, then we can evaluate the performance of a simulation algorithm when it simulates a system that remains in a steady state condition for the entire length of the simulation.
Here we use the fully connected model starting from two different steady state conditions to quantify simulation runtime of a reaction network composed by only slow or fast reactions.These benchmarks are provided in Table 2, third and fourth columns.In both cases, the model has been generated with  = 20, that is, with 20 variables and 380 reactions.All stochastic reaction constants have been set to 1.In the first benchmark (MCM  = 20), the initial state of all model variables has been set to 20 in order to create a network of slow reactions.In the second one (MCM  = 2000), instead, the initial state of all model variables has been set to 2000 in order to create a network of fast reactions.In both conditions, simulation algorithms implemented in HSimulator were the fastest (simulation length 150 time units).
where   and   are random species such that   ∈ { 1 , . . .,   fast } and   ∈ { 1 , . . .,   slow }.The stochastic constant rates of fast reactions are some order of magnitude bigger than the rates of the slow reactions to emphasize the difference of speed between the two subnetworks.The fully connected model has been considered to evaluate the performance of simulation algorithms when a massaction biochemical reaction network of reactions that are all slow (MCM  = 20) or fast (MCM  = 2000) is simulated.Conversely, the multiscale model allowed testing the intermediate condition, where the biochemical network can be divided into two subnetworks working at different time scales.
In Table 2 two benchmarks related to this model are provided.The first one (fifth column, MSM (10, 50)) is related to a model with  fast = 10 and  slow = 50 (60 species and 140 reactions), stochastic reaction constants equal to 10.0 (fast reactions) and 0.001 (slow reactions), initial values of fast species set to 2000, and initial values of slow species set to 20.The second benchmark (last column of the table, MSM (20, 100)) uses the same parameters of the first one except for  fast = 20 and  slow = 100 (120 species and 480 reactions).The provided benchmarks show that the gain in simulation runtime provided by HSimulator is scalable with respect to the complexity of the model (simulation length 10 time units).

Conclusions
We presented HSimulator, a state-of-the-art multithread Java simulator for mass-action well-stirred biochemical reaction systems.HSimulator provides a suite of published state-ofthe-art simulation algorithms including, in the same package, the exact algorithm RSSA and the first publicly available implementation of its hybrid version HRSSA.The benchmarks in the paper show that the HSimulator implementation is often faster than the state-of-the-art simulator COPASI [22].This could open new perspectives in computational systems biology, where often scientists have to balance the accuracy of their simulations with the need of considering large reaction networks modeling complex diseases or disorders.

3. 4 .
The Multiscaled Model.The multiscaled model is a theoretical model with a reaction network that can be divided into two subnetworks working at different time scales.It has been specifically considered to test hybrid simulation strategies because it allows testing the intermediate condition (network made of both fast and slow reactions) between the two scenarios previously considered with the fully connected model (network made of only slow or fast reactions).The first subnetwork of the model is given by a fully connected model of  fast species  1 , . . .,   fast modified by a set of fast reactions.The model is then extended with  slow species  1 , . . .,   slow modified by a set of  slow reactions of the form   +   →   ,  = 1, 2, . . .,  slow , 1), with a stochastic reaction constant   associated with each reaction   ; an initial state X( 0 ) defining molecule abundances at time  0 ; the last time instant  end to be simulated.Output: a time series of states X(),  ∈ [ 0 ;  end ], providing the dynamics of the system.

Table 2 :
Simulation running times of HSimulator and of the state-of-the-art simulator COPASI.Except for the deterministic simulation, HSimulator is demonstrated to be faster in all the considered scenarios (FCM indicates the fully connected model; MSM indicates the multiscaled model).The considered models were simulated for 150 time units (MAPK and fully connected model), 12 time units (Gemcitabine), and 10 time units (multiscaled model).