I have been applied Gillespie algorithm for about three years to simulate stochastic chemical reactions, or more coarse-grained gene regulatory networks. I think it is a good point to check Gillespie’s important original papers and write a report about some of my thoughts…

Here, I used the model of excitable gene regulatory network I discussed in my undergraduate thesis as an example to illustrate how is the significance of fluctuations near chemical instability captured by Gillespie Algorithm.

Given that the algorithm is so well-known and basic for biophysics researchers, the computation steps which I have been taught twice in courses will not be included (but links to example codes below ), but only the physics assumptions behind I learnt from Gillespie’s paper by myself is included.

Ver0.0 2023.05.15 Ver1.0 2023.08.07

An exact method is presented for numerically calculating, within the framework of the stochastic formulation of chemical kinetics, the time evolution of any spatially homogeneous mixture of molecular species which interreact through a specified set of coupled chemical reaction channels. The method is a compact, computer-based, Monte Carlo simulation procedure. 1

On Gillespie Algorithm

D.T. Gillespie’s 1976 paper 1 described a mile-stone computer algorithm to simulate the Markov process of stochastic chemical reactions, which is fully equivalent to the master equation. It is very powerful for theorists to check their derived analytical statistical properties by numerical results.

The paper is elegant and enjoyable to read (and easy!!!). Here I used the example in the paper (see below) to quickly summarize the physical ideas (just for my own understanding, for you, it is better to read the original paper ). My aim for this post is to combine it with a 2014 paper. Suppose there is a reaction

Deterministic Approach and its Limitations

For \(N\) molecular species, a system of \(N\) ODEs describes the time evolution of molecular concentration. The reaction constants are viewed as reaction “rates”.

In the example, \(x_i\) represent the concentrations of \(S_i\), and the ODEs are

where \(k\) is the reaction rate and compared with the stochastic formulation, approximations are involved in the deterministic approach (see below).

Physical Basis of Gillespie Algorithm

The only “assumption” is that the reaction parameter \(c_\mu\), which characterizes reaction \(R_\mu\), can be defined as follows.

Figure from 2.

Apart from the species involved (\(S_i\)) and the form of reactions (\(R_\mu\)), we need to specify the reaction parameter (\(c_\mu\)). How \(c_\mu\) are specified was described elegantly in text.

In short, given spatial homogeneity, the average probability that a particular l-2 molecular pair will collide in the next time interval \(\delta t\) is, to first order in \(\delta t\), \(\langle V_{coll}/V\rangle=\pi d_{12}\langle v_{12}\rangle\delta t/V\), where the brackets denote an average over the velocities of all \(S_1-S_2\) molecular pairs. The condition that nonreactive molecular collisions occur much more frequently than reactive molecular collisions is thus a convenient criterion for applicability of the stochastic formulation of chemical kinetics.

Given thermal (not chemical) equilibrium at temperature \(T\), \(\langle v_{12}\rangle\) can be obtained from Maxwell-Boltzmann distributions \(\langle v_{12}\rangle=(8kT/\pi m_{12})^{1/2}\), where \(m_{12}\) is the reduced mass. Consider only when the kinetic energy of collision exceed the activation energy \(u_u^*\), the reaction occurs, we get \(c_\mu=V^{-1}\pi d_{12}^2(8kT/\pi m_{12})^{1/2}exp(-u_u^*/kT)\).

Gillespie Algorithm, Compared with the Deterministic Approach

In the example, represent the abundances of species are \(X_i\), there will be \(X_1X_2\) distinct combinations of reactant molecules inside \(V\). Thus, the average rate of \(R_\mu\) that occurs inside \(V\) is \(\langle X_1X_2\rangle c_\mu\).

The relationship with reaction rate in the deterministic approach is \(k_\mu=\langle x_1x_2\rangle Vc_\mu/\langle x_1\rangle\langle x_2\rangle\).

For “ordinary” chemical systems in which fluctuations and correlations play no significant role, the method stands as an alternative to the traditional procedure of numerically solving the deterministic reaction rate equations. For nonlinear systems near chemical instabilities, where fluctuations and correlations may invalidate the deterministic equations, the method constitutes an efficient way of numerically examining the predictions of the stochastic master equation. 1

In deterministic approach, we take \(k_\mu\doteq Vc_\mu\). It works well when \(\langle x_1x_2\rangle=\langle x_1\rangle\langle x_2\rangle\) is almost true.

The following part of mathematical derivation from probability is so basic in course that I don’t bother to include.

Probability, Markov Chains and Monte Carlo Method

This part was elegantly explained in the lecture note my supervisor gave me (Johan Paulsson, 10/15/2014, Systems Biology 200), and I don’t think I have any better ideas to rewrite. I just summarize some key points. (I don’t know whether the staff I show here violate copyright. If so, I will sincerely すみません…)

Regard the chemical reaction as the random walk on the discrete lattice on the phase space. If we indicate the probability at state \(A\) as \(P(A,t)\), then we have \(P(A,t+dt=P(A,t)-P_{jump}(t,t+dt)\), while \(P_{jump}(t,t+dt)=P_{jump}(t,t+dt\vert A,t)P(A,t)\), then \(\frac{dP(A,t)}{dt}=-\lim\limits_{dt→0}\frac{P_{jump}(t,t+dt\vert A,t)}{dt}P(A,t)\).

From above, we have \(\lim\limits_{dt→0}\frac{P_{jump}(t,t+dt\vert A,t)}{dt}=c_\mu=const\). I try not to proceed…

Output of Gillespie Algorithm

Simulating the Markovian random walk on the grid of the space of species abundances (the phase space), it can give us the time traces of the time evolution of the system (example here 3a). For complex reaction networks, the analytical solution of the stationary probability distribution on the phase space cannot be written down explicitly, the Gillespie Algorithm gives us the result from numeric computation (example here 4).

On Fluctuations near Instability

Nonlinear Rates and Fluctuations

I would like to show an example of an one-variable dimer degradation to show the relationship between the master equation and the deterministic approach. From Gillespie’s paper 1, we set up a toy model, where the molecule is produced with a constant rate and degraded in a dimer form,

\(x \stackrel{\lambda}{\longrightarrow} x+1 \\x \stackrel{\gamma x(x-1)}{\longrightarrow} x-2\)

Note that nonlinear rates prohibit exact moment closure 5. (In linear case, \(\frac{d\langle x\rangle}{dt}\) contains terms no higher than first moment and you can solve it by given conditions. Then, \(\frac{d\langle x^2\rangle}{dt}\) contains terms no higher than the second moment and you can solve it by the solved first moment… But in the nonlinear case, as I have written, \(\frac{d\langle x\rangle}{dt}\) contains the second order term, and can’t be solved analytically.)

Let’s see something more interesting on my undergrad thesis - there are different types of fixed points in a nonlinear networks!!! Full abstract of the thesis in English can be found here.

The Excitable Gene Regulatory Network

Interlinked positive and negative feedback loops are universal in biological systems, for example, the network of Cyclin B and the network of MAPK. And different kinds of behaviors can be achieved, such as mono-stability, bi-stability, oscillation and excitability. A hallmark of excitability is that a transient input is sufficient for triggering a large response. In electrophysiology, when the potential is raised above a threshold, the neurons can transmit signals. Proteins regulated by genetic networks also fits the features of excitability. For example, when E. coli is under stress, it is likely to switch to the competence from vegetative growth.

The fast positive feedback loop drive the system further away from the stable node/focus, and then the slow negative feedback loop drive the system back to the stable node/focus. (My thesis is to determine the parameter and dynamics of 3-node gene regulatory network which can achieve excitable behavior.)

Fixed Point Analysis

Fixed points analysis in the excitable gene regulatory networks appears in 6, a model of B. subtilis competence. Excitable networks are defined when the ODEs has three fixed points, a stable node/focus, a saddle and an unstable focus. Take the network in Figure.1.b and the ODEs in Box.1 (parameter given by Figure.S3.D), I plotted the three fixed points and the trajectories of the deterministic approach given by the ODEs (the MATLAB code attached 3b!). The saddle defines a manifold. If the initial point start from the “node” side of the stable manifold, it will flow to the node. If start from the “unstable focus” side of the stable manifold, it will flow across the unstable focus and then return to the node, and this is an excitement event. (Whether there will be excitement is totally determined by the initial condition. Similar to the bistable case of Lotka-Volterra competition model, where there are two stable nodes and a saddle…)

### Stochastic Model and Reduced Model

In the case that the stable node/focus is close enough to the saddle (chemical instability), taking stochasticity into account, when the trajectory is bouncing around the stable node/focus, it might cross the saddle’s stable manifold and trigger auto-excitement! This auto-excitement can be captured by Gillespie simulations! I will elaborate it by two examples.

Fluctuation Caused Auto-Excitements Forms a Genetic Timer

The excitable network and auto-excitement was also discussed in 8, where the auto-excitement can be regarded as “a genetic timer” . My math understanding of the fixed points in excitable networks mainly came from this short paper.

This is the figure 1 in the paper. A shows a two-node gene regulatory network consists of interlinked positive and negative feedback loop. B is the deterministic model. The filled red circle is the stable focus, the open red circle is the saddle, and the open red diamond is the unstable focus. If the initial state is on the “excitement” side of the stable manifold defined by the saddle. D is from Gillespie simulations. The saddle is very close of the stable focus. When the trajectory is bouncing due to the stochasticity, it might cross the stable manifold and an auto-excitement is triggered.

Note that, if we adjust the parameter set of the network, making the real-part of the eigenvalue of the unstable focus sufficient close to zero, instead of going to the stable focus directly, the trajectory might circle the unstable focus around for several times before returning back! The time spend for each circle is similar and that’s why it is called the genetic timer!

Transcription and Translation Rates Affects the Level of Intrinsic Noise

In 2002, Ozbudak et al. showed that in a stochastic gene regulatory network model, When the transcription rate is high, variability in protein levels is low, but when the transcription rate is lowered and the translation rate is raised, gene expression is far noisier, even at the same mean. (Summarized in 7)

This statement is explicitly shown in my simulations. Refer to page 18-20 in 3c, the corresponding figures and codes in 3a. In my excitable networks, auto-excitement can be triggered with a larger probability with smaller transcription rates and larger translation rates (at the same mean). Visualized in page 19-20, 3c.

References

[1] Gillespie, D. T. (1976). A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. Journal of computational physics, 22(4), 403-434. online

[2] Gillespie, D. T. (1977). Exact stochastic simulation of coupled chemical reactions. The journal of physical chemistry, 81(25), 2340-2361. online

[3] My undergraduate thesis repository.

• [a] Figure 5.5, Figure 5.6

• [b] Figure 3.2

• [c] Brief introduction

[4] Paper repository

[5] Hilfinger, A., Norman, T. M., Vinnicombe, G., & Paulsson, J. (2016). Constraints on fluctuations in sparsely characterized biological systems. Physical review letters, 116(5), 058101. online

[6] Süel, G. M., Garcia-Ojalvo, J., Liberman, L. M., & Elowitz, M. B. (2006). An excitable gene regulatory circuit induces transient cellular differentiation. Nature, 440(7083), 545-550. online

[7] Raj, A., & Van Oudenaarden, A. (2008). Nature, nurture, or chance: stochastic gene expression and its consequences. Cell, 135(2), 216-226. online

[8] Turcotte, M., Garcia-Ojalvo, J., & Süel, G. M. (2008). A genetic timer through noise-induced stabilization of an unstable state. Proceedings of the National Academy of Sciences, 105(41), 15732-15737. online