MCB111: Mathematics in Biology (Fall 2023)

week 10:

Molecular Population Dynamics as a Markov Process

Stochastic simulation of transcriptional bursting

Sampling fast from an exponential distribution

In the Gillespie algorithm, dwell times \(\tau\) follow a exponential distribution.

The probability density function (PDF) is

\[PDF(t) = W_R e^{-W_R\ t}.\]

In order to sample a dwell time \(\tau\), we used the cumulative probability function (CDF)

\[CDF(\tau) = \int_0^{\tau} W_R\ e^{-W_R\ t}\ dt\]

such that for a random number \(0\leq x\leq 1\), the sampled dwell \(\tau\) time is given by

\[x = CDF(\tau).\]

The CDF of the exponential distribution has a closed form given by

\[CDF(\tau) = \int_0^{\tau} W_R\ e^{-W_R\ t} = 1 - e^{-W_R\ \tau}.\]

Then, the sampled time is given by

\[\tau = -\frac{1}{W_R}\ \log(1-x).\]

If \(x\) is a random number in \([0,1]\), then \(1-x\) is also a random number in the same range, then we can take the sampled dwell time \(\tau\) given by

\[\tau = -\frac{1}{W_R}\ \log(x).\]

Implementing the Gillespie algorithm

For the stochastic process described in Figure 1, the reactions are given by

Figure 1. Regulation of a single gene, transcription and translation, as well as RNA and Protein degradation processes.
r reaction propensity \(w_r\) \(\eta(R)\) \(\eta(P)\)
1 \(\mbox{Gene(I)} \overset{k_b}{\longrightarrow} \mbox{Gene(A)}\) \(k_b\) \(0\) \(0\)
2 \(\mbox{Gene(A)} \overset{k_u}{\longrightarrow} \mbox{Gene(I)}\) \(k_u\) \(0\) \(0\)
3 \(\mbox{Gene(A)} \overset{k_1}{\longrightarrow} RNA\) \(k_1\) \(+1\) \(0\)
4 \(RNA \overset{k_2}{\longrightarrow} \emptyset\) \(k_2 R\) \(-1\) \(0\)
5 \(RNA \overset{k_3}{\longrightarrow} Protein\) \(k_3 R\) \(0\) \(+1\)
6 \(Protein \overset{k_4}{\longrightarrow} \emptyset\) \(k_4 P\) \(0\) \(-1\)

Then the Gillespie algorithm works as follows

\[t\rightarrow t+\tau.\]