Cell Cycle Modeling for Budding Yeast with Stochastic Simulation Algorithms

Background: For biochemical systems, where some chemical species are represented by small numbers of molecules, discrete and stochastic approaches are more appropriate than continuous and deterministic approaches. The continuous deterministic approach using ordinary differential equations is adequate for understanding the average behavior of cells, while the discrete stochastic approach accurately captures noisy events in the growth-division cycle. Since the emergence of the stochastic simulation algorithm (SSA) by Gillespie, alternative algorithms have been developed whose goal is to improve the computational efficiency of the SSA. Results: This paper explains and empirically compares the performance of some of these SSA alternatives on a realistic model. The budding yeast cell cycle provides an excellent example of the need for modeling stochastic effects in mathematical modeling of biochemical reactions. This paper presents a stochastic approximation of the cell cycle for budding yeast using Gillespie’s stochastic simulation algorithm. To compare the stochastic results with the average behavior, the simulation must be run thousands of times. Conclusions: Many of the proposed techniques to accelerate the SSA are not effective on the budding yeast problem, because of the scale of the problem or because underlying assumptions are not satisfied. A load balancing algorithm improved overall performance on a parallel supercomputer. Background The cell-division cycle is the sequence of events that take place in a eukaryotic cell leading to its replication. A growing cell replicates all its components and divides them into two daughter cells, so that each daughter has the information and machinery necessary to repeat the process [1]. Mathematical modeling and computational methods are needed to understand complex yeast control 1 systems. Deterministic mathematical modeling for the budding yeast cell cycle gives the average behavior of populations of dividing cells [2]. However, some major regulatory proteins occur in small numbers such that minor changes in timing and reaction rates can have major inputs on outcomes. Thus, the stochastic approach provides more accurate results than does the deterministic one. In addition, when cell cycle controls are compromised by mutation, random fluctuations are important for modeling the effects of the mutants. Therefore, it is desirable to translate a deterministic cell cycle model into a stochastic model, and simulate the model with an appropriate stochastic method. Gillespie’s stochastic simulation algorithm (SSA) ([3], [4]) is a well-known algorithm using Monte Carlo methods to simulate the chemical reactions. The SSA is an asymptotically exact stochastic method to simulate chemical systems, but the SSA is often slow because it simulates every reaction. Since the SSA emerged, there have been many attempts to improve the computational efficiency ([5], [6], [7]), however, the core principles remain the same. One notable attempt to improve the SSA is the tau-leaping method [6]. Tau-leaping attempts to achieve increased computational efficiency by leaping over many fast reactions. The implicit tau-leaping method compensates for difficulty with stiff systems [8]. Stiff systems are characterized by well separated fast and slow time scales in a dynamic system, the fastest of which is stable. Some approaches try to reduce time consumption with different assumptions such as quasi steady-state approximation (QSSA) [9] and total quasi steady-state approximation (tQSSA) [10] for stiff systems. This paper compares computational efficiency and exactness between SSA, tau-leaping, implicit tau-leaping, QSSA, and tQSSA based on numerical experiments first with a model using simple chemical reactions and stiff systems, and then with the budding yeast model. StochKit [11] is used to do stochastic simulation of the budding yeast model. Because StochKit supports various approximate simulation methods based on the SSA such as explicit and implicit tau-leaping methods, the computational efficiency of the approximation methods can be compared easily by using StochKit. Stochastic methods require that the model be cast in terms of population because they consider reactions with individual molecules. The problem is, however, that ODE models are usually based on concentration values. Therefore, the concentration-based model has to be changed into a population-based model to simulate using a stochastic method. Previous work [12] explained the conversion process using JigCell [13] in detail. StochKit [11] is used to do stochastic simulation of the converted budding yeast model. StochKit supports various approximate simulation methods based on the SSA, but only the exact SSA is used to get precise results. Because the SSA simulates every time step, the SSA is much slower than a deterministic simulation. Moreover, the simulation must be run thousands of times to generate enough data to determine the correct distribution of the behavior. Therefore, it is desirable to run many independent SSA simulations in parallel. StochKit supports MPI for parallel SSA runs, but the user assigns jobs for each processor. Sometimes, the 2 processor times for individual runs are quite different. The result is poor parallel efficiency. This paper presents a dynamic load balancing algorithm that improves the parallel efficiency. The SSA and approximation methods are explained in the next section. The budding yeast cell cycle model and the dynamic load balancing algorithm are presented next. Finally new biological results and numerical comparisons are given. Stochastic Simulation Algorithms SSA Suppose a biochemical system or pathway involves N molecular species {S1, ..., SN}. Xi(t) denotes the number of molecules of species Si at time t. People would like to study the evolution of the state vector X(t) = (X1(t), ..., XN (t)) given that the system was initially in the state vector X(t0). Suppose the system is composed of M reaction channels {R1, ..., RM}. In a constant volume Ω, assume that the system is well-stirred and in thermal equilibrium at some constant temperature. There are two important quantities in reaction channels Rj : the state change vector vj = (v1j , ..., vNj), and propensity function aj . vij is defined as the change in the Si molecules’ population caused by one Rj reaction, and aj(x)dt gives the probability that one Rj reaction will occur in the next infinitesimal time interval [t, t+ dt). The SSA simulates every reaction event ([3], [4]). With X(t) = x, p(τ, j|x, t)dτ is defined as the probability that the next reaction in the system will occur in the infinitesimal time interval [t+ τ, t+ τ+dτ), and will be an Rj reaction. By letting a0(x) ≡ ∑M j=1 aj(x), the equation p(τ, j|x, t) = aj(x) exp(−a0(x)τ), can be obtained. On each step of the SSA, two random numbers r1 and r2 are generated from the uniform (0,1) distribution. From probability theory, the time for the next reaction to occur is given by t+ τ , where τ = 1 a0(x) ln( 1 r1 ). The next reaction index j is given by the smallest integer satisfying