Complexity

Volume 2017, Article ID 1232868, 12 pages

https://doi.org/10.1155/2017/1232868

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

^{1}The Microsoft Research–University of Trento Centre for Computational and Systems Biology (COSBI), Piazza Manifattura, No. 1, 38068 Rovereto, Italy^{2}Department of Computer Science, University of Pisa, Pisa, Italy

Correspondence should be addressed to Luca Marchetti; ue.ibsoc@ittehcram

Received 12 May 2017; Revised 18 August 2017; Accepted 19 November 2017; Published 13 December 2017

Academic Editor: Valeri Mladenov

Copyright © 2017 Luca Marchetti et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

#### 1. Introduction

Computational systems biology is becoming a fundamental tool of life-science research, which aims at developing* models* representing biological phenomena and reliable* computational techniques* for their simulation [1–7].

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 are represented in terms of* abundances* (number of molecules) and reactions are defined aswhere the species on the left of the arrow are* reactants* and the ones on the right are* products*. The stoichiometric coefficients and 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* mass-action 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* Rejection-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 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 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–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-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 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.

##### 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 of* reaction probability density function* (pdf):which provides the probability that the next reaction to apply will be in the infinitesimal time interval , given the system state at time .

The* state of the system* at time is represented by the column vectorwhere is the number of molecules of species in the system at time . In the following we will often write for and for to improve the readability of formulas.

The term in (2) is the* propensity* of in the state , while the* total propensity * is the sum of the propensities of all reactions:The reaction propensities are computed from the stochastic reaction rate :where is the number of distinct reactant combinations for in the state [8]. For standard mass-action kinetics, as considered here, it iswhere are the stoichiometric coefficients of according to (1).

Equation (2) shows that the next reaction fires with a discrete probability and its firing time has an exponential distribution with rate . Gillespie introduced the* Direct Method* (DM) for exactly sampling the pdf by applying the inverse transformation:where and are random numbers drawn from a uniform distribution . The pseudocode of the DM algorithm is in Algorithm 1.