Stochastic Simulations in Systems Biology

Introduction

Molecular biology has undergone significant changes in the last decades. With the enormous amount of data generated by the genomic era, the qualitative reasoning about outcomes of experiments showed inherent limitations. Mathematical models of cellular processes, such as signal transduction (Alon, 2007) are becoming an essential aspect of molecular biology. The dynamic interactions between biomolecules are encapsulated in equations of motion, where multiple biochemical parameters determine rates at which biomolecules are synthesized or degraded, at which they associate, dissociate, or are transformed into other biomolecules. Such equations of motion are often systems of ordinary differential equations (ODE) (Hairer et al., 1991), whose state variables are the concentrations of each type of biomolecules – as it is done in classical chemical kinetics (Laidler, 1987). Yet, some cellular components such as transcription factors or mRNA molecules are present in very few copies and the continuous approach of ODEs falls short under these circumstances. Only a discrete approach capturing the stochastic nature of chemical events can properly describe such dynamics (Wilkinson, 2006; McQuarrie, 1967). In this chapter, the formalism and the main algorithms to perform stochastic simulations of such continuous-time Markov jump processes will be explained. Such simulations are concerned with different time and length scales than classical molecular dynamics simulations. The available methods and implementations of the latter are not applicable to capture the reaction dynamics of large ensembles of biomolecules.

The remaining part of the chapter is organized as follows. The next section introduces the traditional law of mass action and the corresponding rate equations. A discussion on the Markov jump process associated with stochastic chemical kinetics follows. Based on that, a Monte Carlo sampling algorithm developed by Gillespie (Gillespie, 2007) along with its many variants is explained. Subsequently, the first approximate algorithm, the τ-leaping method is derived and some algorithmic aspects of it are highlighted. Further assumptions are made in the following section to obtain the chemical Langevin equation and its simulation is briefly discussed. Finally, we draw conclusions and give some outlook regarding this research area in the last section.

The Law of Mass Action

Chemical reactions, and as a consequence most events occurring in cellular processes, can be represented as transformations of molecular species S_i ∈ {S₁, …, S_n}. A reaction ℜ_j, j∈{1,…,m} with a rate constant k_j ∈ ℝ₊, is written in a general form as

(1) where R_ij ∈ ℕ₀ is the stoichiometric coefficient for the species S_i as a reactant and P_ij ∈ ℕ₀ its coefficient as a product (Heinrich and Schuster, 1996).

In the most simplified approach, all reactions of a process follow the law of mass action, originally proposed by Waage and Guldberg (1864): the rate of a reaction is the product of the concentration of the reactants and a constant. Thus the mass action rate function of reaction ℜ_j is given by

(2) with k ≡ (k₁, …, k_m)^T and where the concentration of each species S_i is denoted as z_i(t), i∈{1,…,n}. The state of the system is thus the time function .

With the stoichiometric matrix N ∈ ℤ^n×m, the element of which are N_ij = P_ij – R_ij, that represents the net effect of reaction ℜ_j on the species S_i, the reaction rate equation can be written in vectorial form as

(3)

In case of reactions that do not follow mass action kinetics (i.e. simplifications as Michaelis-Menten (Heinrich and Schuster, 1996)), equation (3) remains valid, but the reaction rate function cannot be expressed as equation (2).