Modeling Stochastic Gene Regulatory Networks Using Direct Solutions of Chemical Master Equation and Rare Event Sampling

1. Introduction

Molecular components in biological systems, such as proteins, DNA, RNAs, and substrates, are often interacting with each other in a complex reaction network to perform certain biological functions. These networks of interacting biomolecules are at the heart of the regulations of many critical cellular processes, from the regulation of gene expression (Arkin, Ross, & McAdams, 1998; Hasty, Pradines, Dolnik, & Collins, 2000; Levin, 2003; McAdams & Arkin, 1997; Ozbudak, Thattai, Kurtser, Grossman, & van Oudenaarden, 2002), signal transduction (Samoilov, Plyasunov, & Arkin, 2005), to the differentiation of stem cells (Ogawa, 1989). The biological networks are intrinsically stochastic due to thermal fluctuations (McCullagh et al., 2009). The intrinsic stochasticity in these cellular processes originates from reactions involving small copy numbers of molecules. It frequently occur in a cell when molecular concentrations are in the range of 0.1 to 1 nM (typically from about 10 to 100 copies in a cell) (Arkin et al., 1998). For example, the regulation of transcriptions depends on the binding of often a few proteins to a promoter site; the synthesis of protein peptides on ribosome involves a small copy number of molecules; and patterns of cell differentiation depend on initial small copy number events. In these biological processes, fluctuations due to the stochastic behavior intrinsic in small copy number events play important roles.

The importance of stochasticity in cellular functions has been well recognized (Mettetal, Muzzey, Pedraza, Ozbudak, & van Oudenaarden, 2006; Ozbudak et al., 2002; Paulsson & Ehrenberg, 2000; Volfson et al., 2006; Zhou, Chen, & Aihara, 2005). Studies of network models show that stochasticity is important for magnifying signal, sharpening discrimination, and inducing multistability (Paulsson & Ehrenberg, 2000). Understanding the stochastic nature and its consequences for cellular processes involving molecular species of small copy numbers in a network is an important problem. A fundamental framework for studying the full stochasticity is the discrete chemical master equation (dCME). Under this framework, the combination of copy numbers of molecular species defines the microscopic state of the molecular interactions in the network. By explicitly treating microscopic states of reactants, reactions can all be effectively modeled as transitions between microstates, with transition rates determined by the physical properties of the molecules and the cell environment. The reaction trajectories can be modeled in the framework of the dCME, which describes continuous time Markov chains (CTMC) equivalent to the Kolmogorov equation. The probability distribution or probability landscape over these microstates and its time-evolving behavior provide a full description of the properties of a stochastic molecular network.

However, it is challenging to solve a dCME that involves a nontrivial number of species. Analytical solutions of the chemical master equation exist only for very simple cases (Vellela & Qian, 2009). A widely used method is to carry out Monte Carlo simulations of the chemical master equation using the Stochastic Simulation Algorithm (SSA), also known as Gillespie’s algorithm (Gillespie, 1977). Although the SSA approach has found wide applications, it is ineffective in simulating rare events, as most of computing time is spent on following high-probability paths (Daigle, Roh, Gillespie, & Petzold, 2011; Roh, Daigle, Gillespie, & Petzold, 2011; Roh, Gillespie, & Petzold, 2010). A recently developed importance path sampling technique, the Adaptively Biased Sequential Importance Sampling (ABSIS) has been shown to significantly improve rare event sampling efficiency (Cao & Liang, 2013). Alternatively, the chemical master equation can also be approximated using the Fokker-Planck equation (FPE) and the chemical Langevin equation (CLE) (Gillespie, 2000).