Homework 6.1: Noise in a switchable promoter (100 pts)¶
This problem was inspired by Munsky, et al. (Science, 2012). We will explore some aspects of stochastic gene expression using both the analytical and computational techniques that we described in class. In this problem we will contrast two models of stochastic gene expression.
The first model consists of a constitutively active promoter where transcription happens with a constant stochastic rate, while the second model consists of a promoter switching stochastically between an ON and an OFF state. Transcription in this model can only occur when the promoter is in the ON state. This could, for example, describe the open and closed chromatin state. mRNA degradation is present in both models. We will not consider protein in this problem.
The models may be written as chemical equations,
\begin{align} \require{mhchem} &\text{Constitutive: } \ce{ON ->[\beta] mRNA ->[\gamma] \varnothing} \\ &\text{Regulated: } \ce{OFF <=>[k_\mathrm{on}][k_\mathrm{off}] ON ->[\beta] mRNA ->[\gamma] \varnothing}. \end{align}
When the promoter is on, transcription proceeds at a rate \(\beta\). Otherwise, no transcripts are made. Transcripts degrade with a rate \(\gamma\).
a) Write down a master equation describing the temporal dynamics of the probability distribution of the number of RNA transcripts, \(P(n,t)\), for the constitutive model. Nondimensionalize time using the decay rate \(\gamma\).
b) We have seen how to start with a master equation and derive an ODE for the dynamics of the mean number of transcripts, \(\langle n(t) \rangle\). For simple unregulated gene expression, we got
\begin{align} \frac{\mathrm{d}\langle n\rangle}{\mathrm{d}t} = \frac{\beta}{\gamma} - \langle n \rangle, \end{align}
where time is dimensionless (\(t = \gamma t_\mathrm{dimensional}\)). We did this by multiplying both sides of the master equation by \(n\) and then summing over all \(n\).
Now, derive an expression for \(\mathrm{d}\langle n^2 \rangle/\mathrm{d}t\). Next, derive an expression for \(\mathrm{d} \sigma^2 / \mathrm{d}t\). This will depend on both \(\sigma^2(t)\) and \(\langle n(t) \rangle\). Hint:
\begin{align} \frac{\mathrm{d}\sigma^2}{\mathrm{d}t} = \frac{\mathrm{d}}{\mathrm{d}t}\left(\langle n^2 \rangle - \langle n \rangle^2\right) = \frac{\mathrm{d}\langle n^2 \rangle}{\mathrm{d}t} - 2\langle n \rangle\, \frac{\mathrm{d}\langle n \rangle}{\mathrm{d}t}. \end{align}
For initial conditions \(n(0) = \sigma^2(0) = 0\) (i.e., no transcripts exist), show that the Fano factor is unity for all time by verifying that
\begin{align} \langle n(t) \rangle = \sigma^2(t) = \frac{\beta}{\gamma}\left(1 - \mathrm{e}^{-t}\right). \end{align}
c) In this part of the problem, you will perform stochastic simulations of the constitutive model. Start by writing down a system of stochastic chemical reactions for the constitutive model of gene expression. Make sure to enumerate all species, all reactions, and the propensity functions and associated parameters. Then, using the Gillespie stochastic simulation algorithm, implement a simulation of the constitutive model of gene expression. Use \(\beta/\gamma = 10\). Initialize your simulation with no copies of mRNA being present.
Plot snapshots of the simulated mRNA distribution at various time points. Choose the final time such that the system has enough time to settle down to its steady state distribution. Comment on what you see.
Compute the mean, variance, Fano factor, and noise \(\eta\), as a function of time from your SSA samples. Do the computed results match what you obtained in part (b)?
d) Now, we will move on to the regulated model with a switchable promoter. Write down a similar master equation for the regulated case. Now, each “state” is not just the mRNA copy number, but also the state of the promoter (on or off). We will define a variable \(a\in\{0,1\}\) describing this, where \(a = 0\) means the promoter is off and \(a=1\) means it is on. Thus, write a master equation for \(P(a,n,t)\). Again, be sure to nondimensionalize time using \(\gamma\).
e) Using the stochastic simulation algorithm, implement a simulation of the regulated model of gene expression. Start with no transcripts with the promoter being off. Again, use \(\beta/\gamma = 10\) and make sure use to use a sufficient number of SSA samples and run the simulation long enough to reach the steady state distribution. You can try this for various values of \(k_\mathrm{on}\) and \(k_\mathrm{off}\).
It can be shown that the steady state, the mean and Fano factor for this model are, respectively,
\begin{align} &\langle n \rangle = \frac{\beta}{\gamma}\,\frac{k_\mathrm{on}}{k_\mathrm{on} + k_\mathrm{off}},\\ &F = \frac{\sigma^2}{\langle n \rangle} = 1 + \frac{\beta}{\gamma}\,\frac{k_\mathrm{off}}{k_\mathrm{on} + k_\mathrm{off}} \,\frac{1}{1 + (k_\mathrm{on}/\gamma) + (k_\mathrm{off}/\gamma)}. \end{align}
(You can derive this result if you like, but it is not required.) Do your steady state results match this? This result may also serve as a guide for interesting values of \(k_\mathrm{off}\) and \(k_\mathrm{off}\) to try in your SSA calculations.
What is the significance that the Fano factor of the regulated model can never be less than unity? Describe physically what is happening when the Fano factor is large.
Can you get qualitatively different steady state distributions of RNA copy number by varying \(k_\mathrm{on}\) and \(k_\mathrm{off}\)? Plot representative distributions.