# HSimulator: Hybrid Stochastic/Deterministic Simulation of Biochemical Reaction Networks

Authors: Luca Marchetti; Rosario Lombardo; Corrado Priami
Journal: Complexity
 (2017)

---

## Abstract

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).


---

## Body

## 1. Introduction

Computational systems biology is becoming a fundamental tool of life-science research, which aims at developingmodels representing biological phenomena and reliablecomputational techniques for their simulation [1–7].We introduceHSimulator, 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 S1,…,Sn are represented in terms ofabundances (number of molecules) and reactions R1,…,Rm are defined as(1)Rj:v1jS1+⋯+vnjSn⟶v1j′S1+⋯+vnj′Sn,cj,where the species on the left of the arrow arereactants and the ones on the right areproducts. The stoichiometric coefficients vij and vij′ indicate how many molecules of reactant are consumed and how many molecules of product are produced, respectively. The constant cj at the end of the reaction is thestochastic reaction constant introduced by Gillespie [8] for computing reaction firing according to the definition ofmass-action propensities.From the milestone work of Gillespie [8], where theDirect 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 theNext Reaction Method (NRM, [9]) and later theRejection-based Stochastic Simulation Algorithm (RSSA, [10]). The latter is tailored for complex biochemical reactions with time-consuming 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 ofapproximate simulation strategies, which sacrifice accuracy to decrease runtime. Such strategies progressively reduce the stochasticity of the dynamics [11, 12] until reaching adeterministic simulation, where the model is translated to a system ofordinary differential equations (ODEs, [7, 13]) and the dynamics provides an averaged behavior.A drawback of an approximate simulation is thatall 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, ahybrid simulation strategy is able to partition reactions into subsets that are simulated by different algorithms [7, 14, 15]. Hybrid simulation therefore needs tosynchronize the progress of different simulation strategies in order to guarantee the exactness of the simulation of slow reactions. This requires consideringtime-varying transition propensities, which in turn require the computation of integrals to calculate the exact firing time of slow reactions [7, 14, 16–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.TheHybrid 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-of-the-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.


## 2. Materials and Methods

HSimulator provides an implementation of five state-of-the-art 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 Section2.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.


### 2.1. 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 ofreaction probability density function (pdf):(2)Pτ,j∣x,t=ajxe-a0xτ,which provides the probability Pτ,j∣x,tdτ that the next reaction to apply will be Rj in the infinitesimal time interval [t+τ;t+τ+dτ], given the system state x at time t.Thestate of the system at time t is represented by the column vector(3)Xt=X1t,…,XntT,where Xi(t) is the number of molecules of species Si in the system at time t. In the following we will often write x for X(t) and xi for Xi(t) to improve the readability of formulas.The termaj(x) in (2) is thepropensity of Rj in the state x, while thetotal propensitya0(x) is the sum of the propensities of all reactions:(4)a0x=∑j=1majx.The reaction propensities aj(x) are computed from the stochastic reaction rate cj:(5)ajx=cj,ifhjx=0hjxcj,otherwise,where hj(x) is the number of distinct reactant combinations for Rj in the state x [8]. For standard mass-action kinetics, as considered here, it is(6)hjx=∏ixivij=∏ixi!vij!xi-vij!,where vij are the stoichiometric coefficients of Rj according to (1).Equation (2) shows that the next reaction Rj fires with a discrete probability aj(x)/a0(x) and its firing time τ has an exponential distribution with rate a0(x). Gillespie introduced theDirect Method (DM) for exactly sampling the pdf Pτ,j∣x,t by applying the inverse transformation:(7)τ=-ln⁡r1a0x,j=the  smallest  index  s.t.∑i=1jaixa0x≥r2,where r1 and r2 are random numbers drawn from a uniform distribution U(0,1). The pseudocode of the DM algorithm is in Algorithm 1.Algorithm 1:Gillespie’s Direct Method (DM).
Input: a reaction system with n species and m reactions as the one in (1), with a stochastic
reaction constantcj associated with each reaction Rj; an initial state X(t0) defining molecule
abundances at timet0; the last time instant tend to be simulated.
Output: a time series of states X(t),t∈[t0;tend], providing the dynamics of the system.
Pseudocode:
(0)  t≔t0;
(1)  while  t<tend  do
(2)  For each reaction Ri,i=1,…,m, compute reaction propensities ai(x) according
to (5);
(3)  Compute a0(x) according to (4);
(4)  Generate two random numbers r1,r2 in U(0,1);
(5)  τ≔-ln⁡(r1)/a0(x);
(6)  Select Rj such that j satisfies
∑i=1j-1ai(x)a0(x)<r2≤∑i=1jai(x)a0(x);
(7)  Compute X(t+τ) by applying Rj to x;
(8)  t≔t+τ;
(9)  end while


### 2.2. 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–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.Algorithm 2:Rejection-based Stochastic Simulation Algorithm (RSSA).
Input: The same as DM (Algorithm 1) together with a simulation parameter 0<δ<1 for
calculating the fluctuation interval of the system state.
Output: a time series of states X(t),t∈[t0;tend], providing the dynamics of the system.
Pseudocode:
(0)  t≔t0;
(1)  while  t<tend  do
(2)  Define the fluctuation interval X=x_;x¯ according to (8);
(3)  For each reaction Ri,i=1,…,m, compute propensity bounds ai(x_), ai(x¯) and
the total propensitya0(x¯);
(4)  while  t<tend∧x∈X  do
(5)  Generate three random numbers r1,r2,r3 in U(0,1);
(6)  τ:=-ln⁡(r1)/a0(x¯);
(7)  Select Rj such that j satisfies
∑i=1j-1ai(x¯)a0(x¯)<r2≤∑i=1jai(x¯)a0(x¯);
(8)  accepted≔false;
(9)  if  r3≤aj(x_)/aj(x¯)  then
(10)  accepted≔true;
(11)  else
(12)  Compute aj(x) according to (5);
(13)  if  r3≤aj(x)/aj(x¯)  then  accepted≔true;
(14)  end  if
(15)  if  accepted  then
(16)  Compute X(t+τ) by applying Rj to x;
(17)  end  if
(18)  t≔t+τ;
(19)  end  while
(20)  end  whileRSSA computes for each reaction a lower bound and an upper bound of the propensity, which encompass all possible values of reaction propensities over the fluctuation interval. Such bounds are derived by defining afluctuation intervalX=x_;x¯ for the system state x (step (2)), where(8)x_=x-δx,x¯=x+δxand 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 Ri, for i=1,…,m, correspond to the interval [ai(x_),ai(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 reactionRj with probability aj(x¯)/a0(x¯) (step (7)), where a0(x¯)=∑j=1m‍aj(x¯). The candidate reaction is then validated according to a rejection-based procedure that implements the toss of a biased coin with success probability aj(x)/aj(x¯) (steps (8)–(14)). To do this, the algorithm generates a random number r in U(0,1) and moves in one of the following three cases:(1)
Ifr≤aj(x_)/aj(x¯), then Rj is accepted without computing its propensity aj(x) because aj(x_)/aj(x¯)≤aj(x)/aj(x¯) (quick accept);(2)
Ifr>aj(x_)/aj(x¯), then reaction propensity aj(x) is computed and Rj is accepted if r≤aj(x)/aj(x¯) (slow accept);(3)
Ifr>aj(x)/aj(x¯), then Rj is rejected (rejection).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 statex 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.


### 2.3. The Forward Euler Algorithm

The forward Euler algorithm is the simplest example ofdeterministic simulation algorithm. Here the mass-action model is translated to a system ofordinary differential equations (ODEs, [7, 13]) and the dynamics provides an averaged behavior of the system.Consider a biochemical reaction system withS1,…,SN species, R1,…,RM reactions defined according to (1). The mass-action deterministic rate constant kj of Rj can be obtained from the stochastic one cj as(9)kj=cjNAVOrderj-1∏i=1Nvij!,where NA indicates Avogadro’s number, V is the volume where the reaction is taking place, and Orderj is the overall order of reaction Rj, 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 Si, i=1,…,N, is(10)dSidt=∑j=1Mkjvij′-vij∏l=1NSlvlj,i=1,…,N,where squared brackets are used to indicate molar concentration of species ([Si]=Si/(NAV),i=1,…,N). Once the model has been translated to a system of ODEs(11)dXdt=Ft,X,a first-order approximation of the next system state can be computed by(12)Xn+1=Xn+h·Ftn,Xn,where h is a user-defined discretization stepsize. The pseudocode of the Forward Euler algorithm is provided in Algorithm 3.Algorithm 3:Forward Euler method.
Input: a system of ODEs d[X]/dt=F(t,[X]) corresponding to a biochemical reaction
system, the initial state[X(t0)] of the system with species concentrations at time 0, the last
time instanttend to be simulated and the discretization stepsize h.
Output: a time series of states [X(t)],t∈[t0;tend], providing the dynamics of the system in
terms of molar concentrations with discretization stepsizeh.
Pseudocode:
(0)  initialize time t=0 and state [X]=[X(t0)];
(1)  while  t<tend  do
(2)  update [X]=[X]+h⋅F(t,[X]);
(3)  update t=t+h;
(4)  end while


### 2.4. Runge-Kutta-Fehlberg Algorithm (RK45)

The forward Euler algorithm introduced in the previous section is a first-order 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 theRunge-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(13)Xn+1=Xn+25216K1+14082565K3+21974104K4-15K5coupled with a fifth-order RK method(14)X~n+1=Xn+16135K1+665612825K3+2856156430K4-950K5+255K6.Given a system of ODEs as in (11) and a discretization stepsize h, the values of K1,…,K6 are shared between (13) and (14) and they can be computed as(15)K1=h·Ftn,XnK2=h·Ftn+h4,Xn+14K1K3=h·Ftn+3h8,Xn+332K1+932K2K4=h·Ftn+12h13,Xn+19322197K1-72002197K2+72962197K3K5=h·Ftn+h,Xn+439216K1-8K2+3680513K3-8454104K4K6=h·Ftn+h2,Xn-827K1+2K2-35442565K3+18594104K4-1140K5.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:(16)Δn+1=max⁡X~n+1-Xn+1h,where max provides the maximum value along the vector of truncation errors. The error estimate Δn+1 is then compared to the error threshold ϵt specified by the user. When Δn+1≤ϵt, the local truncation error is assumed to be smaller than the threshold, the state [Xn+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 h. In both cases, the value of h is updated as(17)hn+1=hnσ(18)σ=ϵt2Δn+11/4≈0.84ϵtΔn+11/4.When the estimations [Xn+1] and [X~n+1] agree to more significant digits than required, then σ becomes greater than 1 and h is increased. Equation (18) is derived from the general formula σ=(ϵt/Δn+1)1/p, which defines how to update the value of h of an adaptive one-step numerical method of order p, by considering the error estimate Δn+1 and the user-defined threshold ϵt. The additional multiplicative factor of 0.84 is an empirical number commonly added in RK45 implementation to reduce the variability of h, because very high values of the stepsize increase the probability of repeating the next computed step.The implementation of the algorithm is in Algorithm4. The value of h is updated at each step starting from an user-provided initial value h0. The next computed state of the system is accepted only when the estimate of the local truncation error Δ remains below the user-provided threshold ϵt. Steps (16)–(21) are additional steps added to the implementation in order to avoid very large modifications of h in a single step.Algorithm 4:The Runge-Kutta Fehlberg (RK45) algorithm.
Input: a system of ODEs d[X]/dt=F(t,[X]) corresponding to a biochemical reaction
system, the initial state[X(t0)] of the system with species concentrations at time 0, the last
time instanttend to be simulated, an initial value for the discretization stepsize h0 and a
threshold for the maximum local truncation errorϵt.
Output: a time series of states [X(t)],t∈[t0;tend], providing the dynamics of the system in
terms of molar concentrations.
Pseudocode:
(0)  initialize time t=0, state [X]=[X(t0)] and discretization stepsize h=h0;
(1)  while  t<tend  do
(2)  compute K1=h⋅F(t,[X]);
(3)  compute K2=h⋅Ft+h/4,[X]+1/4K1;
(4)  compute K3=h⋅Ft+3h/8,[X]+3/32K1+9/32K2;
(5)  compute K4=h⋅Ft+12h/13,[X]+1932/2197K1-7200/2197K2+7296/2197K3;
(6)  compute K5=h⋅Ft+h,[X]+439/216K1-8K2+3680/513K3-845/4104K4;
(7)  compute K6=h⋅Ft+h/2,[X]-8/27K1+2K2-3544/2565K3+1859/4104K4-11/40K5;
(8)  compute [Xnew]=[X]+25/216K1+1408/2565K3+2197/4104K4-1/5K5;
(9)  compute [X~new]=[X]+16/135K1+6656/12825K3+28561/56430K4-9/50K5+2/55K6;
(10)  compute Δ=max⁡[X~new]-[Xnew]/h according to (16);
(11)  if  Δ≤ϵt  then
(12)  update [X]=[Xnew];
(13)  update t=t+h;
(14)  end if
(15)  compute σ=0.84ϵt/Δ1/4;
(16)  if  σ<0.1  then
(17)  update σ=0.1;
(18)  end if
(19)  if  σ>4  then
(20)  update σ=4;
(21)  end if
(22)  update h=h⋅σ;
(23)  end whileEven 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 ofh often dramatically decreases the simulation runtime.


### 2.5. Hybrid Rejection-Based Stochastic Simulation Algorithm (HRSSA)

A drawback of approximate simulation is thatall 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, ahybrid 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 thesynchronization between the employed simulation strategies in order to guarantee the 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 aj of a slow reaction Rj∈Rslow 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 considertime-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–20], the pdf of the next firing of a slow reaction Rj∈Rslow finally becomes(19)Pslowτ,j∣x,t=ajXt+τe-∫tt+τa0slowXt′dt′,where(20)a0slowx=∑Rj∈Rslowajx.The firing timeτ of the next slow reaction Rj is thus obtained by solving the equation(21)∫tt+τa0slowXt′dt′=-ln⁡r,where r is a random number from U(0,1). Solving (21) is computationally challenging because the state X(t) is changed by fast reactions during time interval [t,t+τ]. 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 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.Algorithm 5:The dynamic reaction partitioning of HRSSA (please refer to [21] for details).
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 Rslow and Rfast of slow and fast reactions, respectively.
Pseudocode:
(0)  Rslow≔∅; Rfast≔∅;
(1)  for each  reaction Rj
(2)  if   (aj(x_)τfast<θ)  then
(3)  Put reaction Rj in the set Rslow;
(4)  else if   (∃Si, modified by Rj with stoichiometric coefficient vij, such that
x_i<γ⋅vij)  then
(5)  Put reaction Rj in the set Rslow;
(6)  else
(7)  Put reaction Rj in the set Rfast;
(8)  end if
(9)  end forAfter reaction partitioning in the two sets of slow and fast reactions, HRSSA computes the sum of upper propensity boundsa0slow(x¯) of slow reactions as defined in (20). The firing time of a candidate slow reaction is then computed as (22)τ=-ln⁡ra0slowx¯,where r is a random number in U(0,1).Under the hypothesis that the system state will remain inside its fluctuation interval in[t,t+τ], we can consider a0slow(x¯) not dependent on time over [t,t+τ] 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[t,t+τ], 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.Algorithm 6:The hybrid rejection-based stochastic simulation algorithm (HRSSA, please refer to [21] for details). In HSimulator, step (9) is implemented by a deterministic Runge-Kutta numerical method.
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(t),t∈[t0;tend], providing the dynamics of the system.
Pseudocode:
(0)  t≔t0;
(1)  while  t<tend  do
(2)  Define the fluctuation interval X=x_;x¯ according to (8);
(3)  For each reaction Ri,i=1,…,m, compute propensity bounds ai(x_) and ai(x¯);
(4)  Update reaction partitioning (sets Rslow and Rfast) by applying the algorithm
in Algorithm5 according to input parameters γ, θ and τfast;
(5)  Compute a0slow(x¯) according to (20);
(6)  updateNeeded≔false;
(7)  while  t<tend∧¬updateNeeded  do
(8)  τ≔-ln⁡(r)/a0slow(x¯), where r is a random number in U(0,1);
(9)  Compute X(t+τ′) by simulating fast reactions (Rfast), at time steps of
maximum lengthτfast, according to an approximate algorithm (either
stochastic or deterministic), whereτ′ is a value τ′≤τ such that X(t+τ′)∈X;
(10)  if  (τ′=τ)  then
(11)  Select a candidate slow reaction Rj∈Rslow by applying RSSA steps (7)–(14)
in Algorithm2;
(12)  if  Rj is accepted by RSSA  then  update X(t+τ′) computed at step (9) by
applyingRj;
(13)  if  X(t+τ′)∉X  then  updateNeeded≔true;
(14)  else
(15)  updateNeeded≔true;
(16)  end if
(17)  t≔t+τ′;
(18)  end while
(19)  end while


## 2.1. 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 ofreaction probability density function (pdf):(2)Pτ,j∣x,t=ajxe-a0xτ,which provides the probability Pτ,j∣x,tdτ that the next reaction to apply will be Rj in the infinitesimal time interval [t+τ;t+τ+dτ], given the system state x at time t.Thestate of the system at time t is represented by the column vector(3)Xt=X1t,…,XntT,where Xi(t) is the number of molecules of species Si in the system at time t. In the following we will often write x for X(t) and xi for Xi(t) to improve the readability of formulas.The termaj(x) in (2) is thepropensity of Rj in the state x, while thetotal propensitya0(x) is the sum of the propensities of all reactions:(4)a0x=∑j=1majx.The reaction propensities aj(x) are computed from the stochastic reaction rate cj:(5)ajx=cj,ifhjx=0hjxcj,otherwise,where hj(x) is the number of distinct reactant combinations for Rj in the state x [8]. For standard mass-action kinetics, as considered here, it is(6)hjx=∏ixivij=∏ixi!vij!xi-vij!,where vij are the stoichiometric coefficients of Rj according to (1).Equation (2) shows that the next reaction Rj fires with a discrete probability aj(x)/a0(x) and its firing time τ has an exponential distribution with rate a0(x). Gillespie introduced theDirect Method (DM) for exactly sampling the pdf Pτ,j∣x,t by applying the inverse transformation:(7)τ=-ln⁡r1a0x,j=the  smallest  index  s.t.∑i=1jaixa0x≥r2,where r1 and r2 are random numbers drawn from a uniform distribution U(0,1). The pseudocode of the DM algorithm is in Algorithm 1.Algorithm 1:Gillespie’s Direct Method (DM).
Input: a reaction system with n species and m reactions as the one in (1), with a stochastic
reaction constantcj associated with each reaction Rj; an initial state X(t0) defining molecule
abundances at timet0; the last time instant tend to be simulated.
Output: a time series of states X(t),t∈[t0;tend], providing the dynamics of the system.
Pseudocode:
(0)  t≔t0;
(1)  while  t<tend  do
(2)  For each reaction Ri,i=1,…,m, compute reaction propensities ai(x) according
to (5);
(3)  Compute a0(x) according to (4);
(4)  Generate two random numbers r1,r2 in U(0,1);
(5)  τ≔-ln⁡(r1)/a0(x);
(6)  Select Rj such that j satisfies
∑i=1j-1ai(x)a0(x)<r2≤∑i=1jai(x)a0(x);
(7)  Compute X(t+τ) by applying Rj to x;
(8)  t≔t+τ;
(9)  end while


## 2.2. 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–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.Algorithm 2:Rejection-based Stochastic Simulation Algorithm (RSSA).
Input: The same as DM (Algorithm 1) together with a simulation parameter 0<δ<1 for
calculating the fluctuation interval of the system state.
Output: a time series of states X(t),t∈[t0;tend], providing the dynamics of the system.
Pseudocode:
(0)  t≔t0;
(1)  while  t<tend  do
(2)  Define the fluctuation interval X=x_;x¯ according to (8);
(3)  For each reaction Ri,i=1,…,m, compute propensity bounds ai(x_), ai(x¯) and
the total propensitya0(x¯);
(4)  while  t<tend∧x∈X  do
(5)  Generate three random numbers r1,r2,r3 in U(0,1);
(6)  τ:=-ln⁡(r1)/a0(x¯);
(7)  Select Rj such that j satisfies
∑i=1j-1ai(x¯)a0(x¯)<r2≤∑i=1jai(x¯)a0(x¯);
(8)  accepted≔false;
(9)  if  r3≤aj(x_)/aj(x¯)  then
(10)  accepted≔true;
(11)  else
(12)  Compute aj(x) according to (5);
(13)  if  r3≤aj(x)/aj(x¯)  then  accepted≔true;
(14)  end  if
(15)  if  accepted  then
(16)  Compute X(t+τ) by applying Rj to x;
(17)  end  if
(18)  t≔t+τ;
(19)  end  while
(20)  end  whileRSSA computes for each reaction a lower bound and an upper bound of the propensity, which encompass all possible values of reaction propensities over the fluctuation interval. Such bounds are derived by defining afluctuation intervalX=x_;x¯ for the system state x (step (2)), where(8)x_=x-δx,x¯=x+δxand 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 Ri, for i=1,…,m, correspond to the interval [ai(x_),ai(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 reactionRj with probability aj(x¯)/a0(x¯) (step (7)), where a0(x¯)=∑j=1m‍aj(x¯). The candidate reaction is then validated according to a rejection-based procedure that implements the toss of a biased coin with success probability aj(x)/aj(x¯) (steps (8)–(14)). To do this, the algorithm generates a random number r in U(0,1) and moves in one of the following three cases:(1)
Ifr≤aj(x_)/aj(x¯), then Rj is accepted without computing its propensity aj(x) because aj(x_)/aj(x¯)≤aj(x)/aj(x¯) (quick accept);(2)
Ifr>aj(x_)/aj(x¯), then reaction propensity aj(x) is computed and Rj is accepted if r≤aj(x)/aj(x¯) (slow accept);(3)
Ifr>aj(x)/aj(x¯), then Rj is rejected (rejection).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 statex 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.


## 2.3. The Forward Euler Algorithm

The forward Euler algorithm is the simplest example ofdeterministic simulation algorithm. Here the mass-action model is translated to a system ofordinary differential equations (ODEs, [7, 13]) and the dynamics provides an averaged behavior of the system.Consider a biochemical reaction system withS1,…,SN species, R1,…,RM reactions defined according to (1). The mass-action deterministic rate constant kj of Rj can be obtained from the stochastic one cj as(9)kj=cjNAVOrderj-1∏i=1Nvij!,where NA indicates Avogadro’s number, V is the volume where the reaction is taking place, and Orderj is the overall order of reaction Rj, 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 Si, i=1,…,N, is(10)dSidt=∑j=1Mkjvij′-vij∏l=1NSlvlj,i=1,…,N,where squared brackets are used to indicate molar concentration of species ([Si]=Si/(NAV),i=1,…,N). Once the model has been translated to a system of ODEs(11)dXdt=Ft,X,a first-order approximation of the next system state can be computed by(12)Xn+1=Xn+h·Ftn,Xn,where h is a user-defined discretization stepsize. The pseudocode of the Forward Euler algorithm is provided in Algorithm 3.Algorithm 3:Forward Euler method.
Input: a system of ODEs d[X]/dt=F(t,[X]) corresponding to a biochemical reaction
system, the initial state[X(t0)] of the system with species concentrations at time 0, the last
time instanttend to be simulated and the discretization stepsize h.
Output: a time series of states [X(t)],t∈[t0;tend], providing the dynamics of the system in
terms of molar concentrations with discretization stepsizeh.
Pseudocode:
(0)  initialize time t=0 and state [X]=[X(t0)];
(1)  while  t<tend  do
(2)  update [X]=[X]+h⋅F(t,[X]);
(3)  update t=t+h;
(4)  end while


## 2.4. Runge-Kutta-Fehlberg Algorithm (RK45)

The forward Euler algorithm introduced in the previous section is a first-order 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 theRunge-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(13)Xn+1=Xn+25216K1+14082565K3+21974104K4-15K5coupled with a fifth-order RK method(14)X~n+1=Xn+16135K1+665612825K3+2856156430K4-950K5+255K6.Given a system of ODEs as in (11) and a discretization stepsize h, the values of K1,…,K6 are shared between (13) and (14) and they can be computed as(15)K1=h·Ftn,XnK2=h·Ftn+h4,Xn+14K1K3=h·Ftn+3h8,Xn+332K1+932K2K4=h·Ftn+12h13,Xn+19322197K1-72002197K2+72962197K3K5=h·Ftn+h,Xn+439216K1-8K2+3680513K3-8454104K4K6=h·Ftn+h2,Xn-827K1+2K2-35442565K3+18594104K4-1140K5.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:(16)Δn+1=max⁡X~n+1-Xn+1h,where max provides the maximum value along the vector of truncation errors. The error estimate Δn+1 is then compared to the error threshold ϵt specified by the user. When Δn+1≤ϵt, the local truncation error is assumed to be smaller than the threshold, the state [Xn+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 h. In both cases, the value of h is updated as(17)hn+1=hnσ(18)σ=ϵt2Δn+11/4≈0.84ϵtΔn+11/4.When the estimations [Xn+1] and [X~n+1] agree to more significant digits than required, then σ becomes greater than 1 and h is increased. Equation (18) is derived from the general formula σ=(ϵt/Δn+1)1/p, which defines how to update the value of h of an adaptive one-step numerical method of order p, by considering the error estimate Δn+1 and the user-defined threshold ϵt. The additional multiplicative factor of 0.84 is an empirical number commonly added in RK45 implementation to reduce the variability of h, because very high values of the stepsize increase the probability of repeating the next computed step.The implementation of the algorithm is in Algorithm4. The value of h is updated at each step starting from an user-provided initial value h0. The next computed state of the system is accepted only when the estimate of the local truncation error Δ remains below the user-provided threshold ϵt. Steps (16)–(21) are additional steps added to the implementation in order to avoid very large modifications of h in a single step.Algorithm 4:The Runge-Kutta Fehlberg (RK45) algorithm.
Input: a system of ODEs d[X]/dt=F(t,[X]) corresponding to a biochemical reaction
system, the initial state[X(t0)] of the system with species concentrations at time 0, the last
time instanttend to be simulated, an initial value for the discretization stepsize h0 and a
threshold for the maximum local truncation errorϵt.
Output: a time series of states [X(t)],t∈[t0;tend], providing the dynamics of the system in
terms of molar concentrations.
Pseudocode:
(0)  initialize time t=0, state [X]=[X(t0)] and discretization stepsize h=h0;
(1)  while  t<tend  do
(2)  compute K1=h⋅F(t,[X]);
(3)  compute K2=h⋅Ft+h/4,[X]+1/4K1;
(4)  compute K3=h⋅Ft+3h/8,[X]+3/32K1+9/32K2;
(5)  compute K4=h⋅Ft+12h/13,[X]+1932/2197K1-7200/2197K2+7296/2197K3;
(6)  compute K5=h⋅Ft+h,[X]+439/216K1-8K2+3680/513K3-845/4104K4;
(7)  compute K6=h⋅Ft+h/2,[X]-8/27K1+2K2-3544/2565K3+1859/4104K4-11/40K5;
(8)  compute [Xnew]=[X]+25/216K1+1408/2565K3+2197/4104K4-1/5K5;
(9)  compute [X~new]=[X]+16/135K1+6656/12825K3+28561/56430K4-9/50K5+2/55K6;
(10)  compute Δ=max⁡[X~new]-[Xnew]/h according to (16);
(11)  if  Δ≤ϵt  then
(12)  update [X]=[Xnew];
(13)  update t=t+h;
(14)  end if
(15)  compute σ=0.84ϵt/Δ1/4;
(16)  if  σ<0.1  then
(17)  update σ=0.1;
(18)  end if
(19)  if  σ>4  then
(20)  update σ=4;
(21)  end if
(22)  update h=h⋅σ;
(23)  end whileEven 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 ofh often dramatically decreases the simulation runtime.


## 2.5. Hybrid Rejection-Based Stochastic Simulation Algorithm (HRSSA)

A drawback of approximate simulation is thatall 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, ahybrid 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 thesynchronization between the employed simulation strategies in order to guarantee the 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 aj of a slow reaction Rj∈Rslow 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 considertime-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–20], the pdf of the next firing of a slow reaction Rj∈Rslow finally becomes(19)Pslowτ,j∣x,t=ajXt+τe-∫tt+τa0slowXt′dt′,where(20)a0slowx=∑Rj∈Rslowajx.The firing timeτ of the next slow reaction Rj is thus obtained by solving the equation(21)∫tt+τa0slowXt′dt′=-ln⁡r,where r is a random number from U(0,1). Solving (21) is computationally challenging because the state X(t) is changed by fast reactions during time interval [t,t+τ]. 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 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.Algorithm 5:The dynamic reaction partitioning of HRSSA (please refer to [21] for details).
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 Rslow and Rfast of slow and fast reactions, respectively.
Pseudocode:
(0)  Rslow≔∅; Rfast≔∅;
(1)  for each  reaction Rj
(2)  if   (aj(x_)τfast<θ)  then
(3)  Put reaction Rj in the set Rslow;
(4)  else if   (∃Si, modified by Rj with stoichiometric coefficient vij, such that
x_i<γ⋅vij)  then
(5)  Put reaction Rj in the set Rslow;
(6)  else
(7)  Put reaction Rj in the set Rfast;
(8)  end if
(9)  end forAfter reaction partitioning in the two sets of slow and fast reactions, HRSSA computes the sum of upper propensity boundsa0slow(x¯) of slow reactions as defined in (20). The firing time of a candidate slow reaction is then computed as (22)τ=-ln⁡ra0slowx¯,where r is a random number in U(0,1).Under the hypothesis that the system state will remain inside its fluctuation interval in[t,t+τ], we can consider a0slow(x¯) not dependent on time over [t,t+τ] 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[t,t+τ], 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.Algorithm 6:The hybrid rejection-based stochastic simulation algorithm (HRSSA, please refer to [21] for details). In HSimulator, step (9) is implemented by a deterministic Runge-Kutta numerical method.
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(t),t∈[t0;tend], providing the dynamics of the system.
Pseudocode:
(0)  t≔t0;
(1)  while  t<tend  do
(2)  Define the fluctuation interval X=x_;x¯ according to (8);
(3)  For each reaction Ri,i=1,…,m, compute propensity bounds ai(x_) and ai(x¯);
(4)  Update reaction partitioning (sets Rslow and Rfast) by applying the algorithm
in Algorithm5 according to input parameters γ, θ and τfast;
(5)  Compute a0slow(x¯) according to (20);
(6)  updateNeeded≔false;
(7)  while  t<tend∧¬updateNeeded  do
(8)  τ≔-ln⁡(r)/a0slow(x¯), where r is a random number in U(0,1);
(9)  Compute X(t+τ′) by simulating fast reactions (Rfast), at time steps of
maximum lengthτfast, according to an approximate algorithm (either
stochastic or deterministic), whereτ′ is a value τ′≤τ such that X(t+τ′)∈X;
(10)  if  (τ′=τ)  then
(11)  Select a candidate slow reaction Rj∈Rslow by applying RSSA steps (7)–(14)
in Algorithm2;
(12)  if  Rj is accepted by RSSA  then  update X(t+τ′) computed at step (9) by
applyingRj;
(13)  if  X(t+τ′)∉X  then  updateNeeded≔true;
(14)  else
(15)  updateNeeded≔true;
(16)  end if
(17)  t≔t+τ′;
(18)  end while
(19)  end while


## 3. 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 mass-action 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 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.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).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.Table 1
Simulation parameters used for the benchmarks of the compared tools.
Algorithm
Software
Parameters
DM
HSimulator
No parameters needed
DM
COPASI
No parameters needed
RSSA
HSimulator
δ
=
0.1
LSODA
COPASI
Relative and absolute thresholds:10-3
RK45
HSimulator
Initial stepsize:10-2
Error threshold:10-3
Hybrid Runge-Kutta
COPASI
Maximum number of reactant molecules for slow reactions: 100
RK4 stepsize:10-2
HRSSA
HSimulator
δ
=
0.1; τfast=10-2; γ=100; θ=10The 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 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 Table2. 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 i=20) or fast (MCM i=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 for details). The provided benchmarks seem to indicate that HSimulator is more scalable than COPASI with respect to the growing complexity of the model. This result is interesting, especially from an implementer’s point of view, showing how a recent Java-based implementation (HSimulator), in the given conditions, may compare favourably with respect to C++ (COPASI).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 for150 time units (MAPK and fully connected model), 12 time units (Gemcitabine), and 10 time units (multiscaled model).
Algorithm
MAPK
Gemcitabine
FCMi=20
FCMi=2000
MSM(10,50)
MSM(20,100)
COPASI - DM
2.39 sec.
0.88 sec.
5.09 sec.
415.04 sec.
28.87 sec.
262.1 sec.
HSimulator - DM
2.23 sec. (−6%)
0.7 sec. (−12%)
1.61 sec. (−68%)
168.75 sec. (−59%)
11.0 sec. (−62%)
120.84 sec. (−54%)
HSimulator - RSSA
2.03 sec. (−15%)
0.66 sec. (−25%)
1.10 sec. (−78%)
77.46 sec. (−81%)
4.25 sec. (−85%)
37.23 sec. (−86%)
COPASI - LSODA
0.06 sec.
0.02 sec.
0.56 sec.
0.08 sec.
0.08 sec.
0.21 sec.
HSimulator - RK45
0.06 sec. (+0%)
0.04 sec. (+100%)
0.001 sec. (−99%)
0.001 sec. (−99%)
0.13 sec. (+63%)
0.52 sec. (+148%)
COPASI - Hyb. RK
0.21 sec
0.54 sec.
26.74 sec.
14.67 sec.
1.93 sec.
35.49 sec.
HSimulator - HRSSA
0.1 sec. (−52%)
0.38 sec. (−30%)
1.10 sec. (−96%)
0.8 sec. (−95%)
2.53 sec. (−31%)
5.89 sec. (−83%)


### 3.1. The MAPK Cascade

TheMitogen-Activated Protein Kinase cascade (MAPK cascade) is one of the most important and intensively studied signaling pathways [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 signalE1 activates MAPKKK by phosphorylation, which in turn activates MAPKK. Once activated, MAPKK activates MAPK. When E1 is added to the system, the output of activated MAPK increases rapidly. By removing the signal E1, the output level of activated MAPK reverts back to zero. Reverse reactions are triggered by the signal E2 and by the MAPK phosphatases KKpase and Kpase. The reaction-based representation is(23)R1:KKK+E1⟶KKKp+E1R2:KKKp+E2⟶KKK+E2R3:KK+KKKp⟶KKp+KKKpR4:KKp+KKpase⟶KK+KKpaseR5:KKp+KKKp⟶KKpp+KKKpR6:KKpp+KKpase⟶KKp+KKpaseR7:K+KKpp⟶Kp+KKppR8:Kp+Kpase⟶K+KpaseR9:Kp+KKpp⟶Kpp+KKppR10:Kpp+Kpase⟶Kp+Kpase,where KKK denotes MAPKKK, KK denotes MAPKK, and K denotes MAPK and tags p and pp indicate phosphorylation and double phosphorylation, respectively.A simulation benchmark of the MAPK cascade is in Table2 (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 E1=2000, E2=2000, K=200000, KK=200000, and KKK=20000. All the other species are initialized with 0 molecules.


### 3.2. The Gemcitabine Model

Gemcitabine is a drug for non-small-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 Box1) is transported from plasma into the cell through the cell membrane (dFdCout→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 dFdCTP. Alternatively, it is further phosphorylated. The Gemcitabine triphosphates dFdCTP and dFdUTP 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.Box 1:Reactions of the Gemcitabine model [23].
R
1
:
dFd
C
out
⟶
dFdC
R
2
:
dFdC
⟶
dFd
C
out
R
3
:
dFd
C
out
⟶
dFdU
R
4
:
dFdU
⟶
dFd
U
out
R
5
:
dFdC
⟶
dFdU
R
6
:
dFdC
+
dCK
⟶
dFd
C
MP
+
dCK
R
7
:
dFd
C
MP
⟶
dFdC
R
8
:
dFd
C
MP
⟶
dFd
C
DP
R
9
:
dFd
C
DP
⟶
dFd
C
MP
R
10
:
dFd
C
DP
⟶
dFd
C
TP
R
11
:
dFd
C
TP
⟶
dFd
C
DP
R
12
:
dFdU
+
dCK
⟶
dFd
U
MP
+
dCK
R
13
:
dFd
U
MP
⟶
dFdU
R
14
:
dFd
U
MP
⟶
dFd
U
DP
R
15
:
dFd
U
DP
⟶
dFd
U
MP
R
16
:
dFd
U
DP
⟶
dFd
U
TP
R
17
:
dFd
U
TP
⟶
dFd
U
DP
R
18
:
dFdC
⟶
dFdU
R
19
:
dFd
C
MP
+
dCMPD
⟶
dFd
U
MP
+
dCMPD
R
20
:
∅
⟶
CDP
R
21
:
CDP
+
RR
⟶
dCDP
+
RR
R
22
:
dCDP
⟶
dCTP
R
23
:
dFd
C
DP
+
RR
⟶
dFd
C
DP
:
RR
R
24
:
dFd
C
DP
:
RR
⟶
dFd
C
DP
+
RR
R
25
:
dCTP
+
dCK
⟶
dCTP
:
dCK
R
26
:
dCTP
:
dCK
⟶
dCTP
+
dCK
R
27
:
dFd
C
TP
+
dCMPD
⟶
dFd
C
TP
:
dCMPD
R
28
:
dFd
C
TP
:
dCMPD
⟶
dFd
C
TP
+
dCMPD
R
29
:
dFd
C
TP
⟶
dFd
C
TP
:
DNA
R
30
:
dFd
U
TP
⟶
dFd
U
TP
:
DNA
R
31
:
dCTP
⟶
dCTP
:
DNAA simulation benchmark of the Gemcitabine model is in Table2 (second column, simulation length 12 time units). In the computed benchmark, stochastic reaction constants are [23]: c1=9.97, c2=2.61e-4, c3=4.72e-6, c4=0.05, c5=0, c6=0.001, c7=0.087, c8=2.37, c9=0.21, c10=2.52, c11=1.45, c12=9.68e-4, c13=5.6e-5, c14=0.078, c15=0.004, c16=0.16, c17=9.05e-5, c18=4.76e-4, c19=4.559e-6, c20=1000, c21=0.05, c22=25.20, c23=1e-5, c24=0.1, c25=1e-5, c26=0.1, c27=1e-7, c28=0.1, c29=0.05, c30=7.37e-4, and c31=0.5. Initially all the variables have been set to 0 except for dFdCout=100000, dCK=1000, RR=1000, dCMPD=1000, and CDP=2000.


### 3.3. The Fully Connected Model

The fully connected model is a theoretical model that consists ofn species S1,S2,…,Sn and a set of reactions that reversibly convert each species into each other of the form (24)Si⟷Sj,i,j=1,…,n,i≠j.The fully connected model is useful to evaluate the scalability of algorithms (variation in runtime with growing values of n). 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 Table2, third and fourth columns. In both cases, the model has been generated with n=20, that is, with 20 variables and 380 reactions. All stochastic reaction constants have been set to 1. In the first benchmark (MCM i=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 i=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).


### 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 ofnfast species S1,…,Snfast modified by a set of fast reactions. The model is then extended with nslow species s1,…,snslow modified by a set of nslow reactions of the form(25)Si+sj⟶sk,j=1,2,…,nslow, where Si and sk are random species such that Si∈{S1,…,Snfast} and sk∈{s1,…,snslow}. 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 mass-action biochemical reaction network of reactions that are all slow (MCMi=20) or fast (MCM i=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 Table2 two benchmarks related to this model are provided. The first one (fifth column, MSM (10,50)) is related to a model with nfast=10 and nslow=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 nfast=20 and nslow=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).


## 3.1. The MAPK Cascade

TheMitogen-Activated Protein Kinase cascade (MAPK cascade) is one of the most important and intensively studied signaling pathways [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 signalE1 activates MAPKKK by phosphorylation, which in turn activates MAPKK. Once activated, MAPKK activates MAPK. When E1 is added to the system, the output of activated MAPK increases rapidly. By removing the signal E1, the output level of activated MAPK reverts back to zero. Reverse reactions are triggered by the signal E2 and by the MAPK phosphatases KKpase and Kpase. The reaction-based representation is(23)R1:KKK+E1⟶KKKp+E1R2:KKKp+E2⟶KKK+E2R3:KK+KKKp⟶KKp+KKKpR4:KKp+KKpase⟶KK+KKpaseR5:KKp+KKKp⟶KKpp+KKKpR6:KKpp+KKpase⟶KKp+KKpaseR7:K+KKpp⟶Kp+KKppR8:Kp+Kpase⟶K+KpaseR9:Kp+KKpp⟶Kpp+KKppR10:Kpp+Kpase⟶Kp+Kpase,where KKK denotes MAPKKK, KK denotes MAPKK, and K denotes MAPK and tags p and pp indicate phosphorylation and double phosphorylation, respectively.A simulation benchmark of the MAPK cascade is in Table2 (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 E1=2000, E2=2000, K=200000, KK=200000, and KKK=20000. All the other species are initialized with 0 molecules.


## 3.2. The Gemcitabine Model

Gemcitabine is a drug for non-small-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 Box1) is transported from plasma into the cell through the cell membrane (dFdCout→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 dFdCTP. Alternatively, it is further phosphorylated. The Gemcitabine triphosphates dFdCTP and dFdUTP 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.Box 1:Reactions of the Gemcitabine model [23].
R
1
:
dFd
C
out
⟶
dFdC
R
2
:
dFdC
⟶
dFd
C
out
R
3
:
dFd
C
out
⟶
dFdU
R
4
:
dFdU
⟶
dFd
U
out
R
5
:
dFdC
⟶
dFdU
R
6
:
dFdC
+
dCK
⟶
dFd
C
MP
+
dCK
R
7
:
dFd
C
MP
⟶
dFdC
R
8
:
dFd
C
MP
⟶
dFd
C
DP
R
9
:
dFd
C
DP
⟶
dFd
C
MP
R
10
:
dFd
C
DP
⟶
dFd
C
TP
R
11
:
dFd
C
TP
⟶
dFd
C
DP
R
12
:
dFdU
+
dCK
⟶
dFd
U
MP
+
dCK
R
13
:
dFd
U
MP
⟶
dFdU
R
14
:
dFd
U
MP
⟶
dFd
U
DP
R
15
:
dFd
U
DP
⟶
dFd
U
MP
R
16
:
dFd
U
DP
⟶
dFd
U
TP
R
17
:
dFd
U
TP
⟶
dFd
U
DP
R
18
:
dFdC
⟶
dFdU
R
19
:
dFd
C
MP
+
dCMPD
⟶
dFd
U
MP
+
dCMPD
R
20
:
∅
⟶
CDP
R
21
:
CDP
+
RR
⟶
dCDP
+
RR
R
22
:
dCDP
⟶
dCTP
R
23
:
dFd
C
DP
+
RR
⟶
dFd
C
DP
:
RR
R
24
:
dFd
C
DP
:
RR
⟶
dFd
C
DP
+
RR
R
25
:
dCTP
+
dCK
⟶
dCTP
:
dCK
R
26
:
dCTP
:
dCK
⟶
dCTP
+
dCK
R
27
:
dFd
C
TP
+
dCMPD
⟶
dFd
C
TP
:
dCMPD
R
28
:
dFd
C
TP
:
dCMPD
⟶
dFd
C
TP
+
dCMPD
R
29
:
dFd
C
TP
⟶
dFd
C
TP
:
DNA
R
30
:
dFd
U
TP
⟶
dFd
U
TP
:
DNA
R
31
:
dCTP
⟶
dCTP
:
DNAA simulation benchmark of the Gemcitabine model is in Table2 (second column, simulation length 12 time units). In the computed benchmark, stochastic reaction constants are [23]: c1=9.97, c2=2.61e-4, c3=4.72e-6, c4=0.05, c5=0, c6=0.001, c7=0.087, c8=2.37, c9=0.21, c10=2.52, c11=1.45, c12=9.68e-4, c13=5.6e-5, c14=0.078, c15=0.004, c16=0.16, c17=9.05e-5, c18=4.76e-4, c19=4.559e-6, c20=1000, c21=0.05, c22=25.20, c23=1e-5, c24=0.1, c25=1e-5, c26=0.1, c27=1e-7, c28=0.1, c29=0.05, c30=7.37e-4, and c31=0.5. Initially all the variables have been set to 0 except for dFdCout=100000, dCK=1000, RR=1000, dCMPD=1000, and CDP=2000.


## 3.3. The Fully Connected Model

The fully connected model is a theoretical model that consists ofn species S1,S2,…,Sn and a set of reactions that reversibly convert each species into each other of the form (24)Si⟷Sj,i,j=1,…,n,i≠j.The fully connected model is useful to evaluate the scalability of algorithms (variation in runtime with growing values of n). 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 Table2, third and fourth columns. In both cases, the model has been generated with n=20, that is, with 20 variables and 380 reactions. All stochastic reaction constants have been set to 1. In the first benchmark (MCM i=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 i=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).


## 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 ofnfast species S1,…,Snfast modified by a set of fast reactions. The model is then extended with nslow species s1,…,snslow modified by a set of nslow reactions of the form(25)Si+sj⟶sk,j=1,2,…,nslow, where Si and sk are random species such that Si∈{S1,…,Snfast} and sk∈{s1,…,snslow}. 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 mass-action biochemical reaction network of reactions that are all slow (MCMi=20) or fast (MCM i=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 Table2 two benchmarks related to this model are provided. The first one (fifth column, MSM (10,50)) is related to a model with nfast=10 and nslow=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 nfast=20 and nslow=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).


## 4. 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-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-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.



---

*Source: 1232868-2017-12-13.xml*