cuTauLeaping: A GPU-Powered Tau-Leaping Stochastic Simulator for Massive Parallel Analyses of Biological Systems

Tau-leaping is a stochastic simulation algorithm that efficiently reconstructs the temporal evolution of biological systems, modeled according to the stochastic formulation of chemical kinetics. The analysis of dynamical properties of these systems in physiological and perturbed conditions usually requires the execution of a large number of simulations, leading to high computational costs. Since each simulation can be executed independently from the others, a massive parallelization of tau-leaping can bring to relevant reductions of the overall running time. The emerging field of General Purpose Graphic Processing Units (GPGPU) provides power-efficient high-performance computing at a relatively low cost. In this work we introduce cuTauLeaping, a stochastic simulator of biological systems that makes use of GPGPU computing to execute multiple parallel tau-leaping simulations, by fully exploiting the Nvidia's Fermi GPU architecture. We show how a considerable computational speedup is achieved on GPU by partitioning the execution of tau-leaping into multiple separated phases, and we describe how to avoid some implementation pitfalls related to the scarcity of memory resources on the GPU streaming multiprocessors. Our results show that cuTauLeaping largely outperforms the CPU-based tau-leaping implementation when the number of parallel simulations increases, with a break-even directly depending on the size of the biological system and on the complexity of its emergent dynamics. In particular, cuTauLeaping is exploited to investigate the probability distribution of bistable states in the Schlögl model, and to carry out a bidimensional parameter sweep analysis to study the oscillatory regimes in the Ras/cAMP/PKA pathway in S. cerevisiae.


Supporting Information -Text S1
Stochastic biological models

Michaelis-Menten kinetics
The Michaelis-Menten (MM) model is a simple example of enzyme kinetics consisting in 3 reactions, which describe the binding of an enzyme E to a substrate S to form an intermediate complex ES, that is eventually converted into a final product P [1]. These reactions, together with the values of the associated stochastic constants, are given in Table 1. The initial molecular amounts used in this work, given as number of molecules, are listed in Table 2.

Prokaryotic auto-regulatory gene network
The prokaryotic gene network (PGN) [2] describes an auto-regulation mechanism of gene expression, whereby a gene (DN A) that codes for a protein (P ) is inhibited by binding to a dimer of the protein itself (DN A:P 2 ). Gene expression is a good example of stochasticity in biological systems: the transcriptional regulators are present in a few copies, so that the binding and release of the regulators can be expressed in probabilistic terms. The reactions describing the molecular interactions occurring in PGN, together with the values of the associated stochastic constants, are given in Table  3. The only initial molecular amount used in this work is DN A=200 molecules; all other molecular species are generated in the system as long as reactions are applied.
Remark : λ denotes the degradation of the reactant The Schlögl system [3,4] is one of the simplest prototypes of chemical systems presenting a bistable dynamical behavior, i.e., the capacity of switching between two different stable steady states in response to some chemical signaling (see, e.g., [5][6][7] and references therein). The Schlögl model consists of 4 chemical reactions and 3 molecular species, listed in Table 4. The initial molecular amounts used in this work are given in Table 5.  In the yeast Saccharomyces cerevisiae, the Ras/cAMP/PKA pathway plays a major role in the regulation of metabolism, stress resistance and cell cycle progression [8,9]. This pathway controls more than 90% of all genes that are regulated by glucose through the activation of the protein kinase A (PKA), that is able to phosphorylate a plethora of downstream proteins. PKA is activated by the binding of the second messenger cyclic-AMP (cAMP), which is synthesized by the adenylate cyclase Cyr1. The activity of Cyr1 is controlled by the monomeric GTPases Ras1 and Ras2, which cycle between a GTP-bound active state and a GDP-bound inactive state. In turn, Ras proteins are positively regulated by protein Cdc25, a Ras-GEF (Guanine Nucleotide Exchange Factor) that stimulates the GDP to GTP exchange, and negatively regulated by proteins Ira1 and Ira2, two Ras-GAP (GTPase Activating Proteins) that stimulate the GTPase activity of Ras proteins. The degradation of cAMP is governed by two phosphodiesterases, Pde1 and Pde2. These two enzymes constitute a major negative feedback in this pathway: the low-affinity phosphodiesterase Pde1 is active under the positive regulation of PKA, while the high-affinity phosphodiesterase Pde2 is active in the basal level regulation of cAMP. The reactions describing the interactions occurring in the Ras/cAMP/PKA pathway, together with the values of the associated stochastic constants, are given in Table 6 (see also [10][11][12] for further details). In particular: • reactions r 1 , . . . , r 10 describe the switch cycle of Ras2 protein between its inactive state (Ras2-GDP) and active state (Ras2-GTP), regulated by the activity of the GEF Cdc25 and of the GAP Ira2; • reactions r 11 , r 12 , r 13 describe the synthesis of cAMP through the activation of the adenylate cyclase Cyr1, mediated by Ras2-GTP; • reactions r 14 , . . . , r 25 describe the activation of PKA, mediated by the reversible binding of cAMP to its two regulatory subunits, and the subsequent dissociation of the PKA tetramer, which releases the two catalytic subunits; • reactions r 26 , . . . , r 33 describe the activity of the two phosphodiesterases Pde1 and Pde2, that carry out the degradation of cAMP. The activation of Pde1 is regulated by the catalytic subunits of PKA, and it represents one of the main negative feedback control exerted by PKA within this pathway; • reactions r 34 , r 35 describe the negative feedback exerted by PKA on Cdc25, whose effect is modeled as a partial inactivation of the GEF activity and a reduction of the active state level of Ras2-GTP.
The initial molecular amounts used in this work are summarized in Table 7.