Homework 5

(c) 2019 Justin Bois and Michael Elowitz. With the exception of pasted graphics, where the source is noted, this work is licensed under a Creative Commons Attribution License CC-BY 4.0. All code contained herein is licensed under an MIT license.

This document was prepared at Caltech with financial support from the Donna and Benjamin M. Rosen Bioengineering Center.

This homework was generated from a Jupyter notebook. You can download the notebook here.


5.1: Delay oscillators, 50 pts

In class, we have seen delay oscillators in which a gene is autoinhibitory with delay. In this problem, we will consider a delay oscillator in which we have autoinhibition, but also enzymatic degradation in addition to the normal autodegradation due to dilution. We will do an analytical and numerical treatment to see what we can learn about the dynamics of delay oscillators. This model was proposed by Mather, et al. (PRL 2009).

a) Let $X$ be the concentration of the protein of interest that undergoes the dynamics described above. We can describe its dynamics with the ODE

\begin{align} \frac{\mathrm{d}X}{\mathrm{d}t} = \beta\,\frac{k^n}{k^n + \left(X(t-\tau)\right)^n} - \delta\,\frac{X}{K + X} - \gamma X. \end{align}

Explain what each term means. Why do we write the delay as we do? What physical process might this be describing?

b) Show that the ODE may be nondimensionalized to read

\begin{align} \frac{\mathrm{d}x}{\mathrm{d}t} = \frac{\beta}{1 + \left(x(t-\tau)\right)^n} - \delta\,\frac{x}{\kappa + x} - x, \end{align}

where we have redefined parameters and variables to be dimensionless.

c) Numerically solve the DDE for $\beta = 1000$, $\delta = 200$, $\kappa = 0.1$, $n = 4$. Perform the calculation for $\tau = 0.1$ and $\tau = 10$. From these numerical calculations, compute the amplitude an period of the oscillations.

d) We will now treat the delay oscillator analytically. As we has seen elsewhere in the class, analytical treatments in certain limits (in which analytical progress is tractable) can reveal insights about the functioning of a circuit. First, we will consider purely switch-like autodegradation ($n \to \infty$). This means that the rate of production of protein is $\beta$ if $x(t-\tau) < 1$ and zero otherwise. This approximates the case of highly ultrasensitive repression. Next, we will consider that the enzymatic degradation is always saturated and progresses at a rate $\delta$. This is the case is $\kappa \ll 1$. So, our approximate ODE is

\begin{align} \frac{\mathrm{d}x}{\mathrm{d}t} \approx \beta[1 - \theta(x(t-\tau) - 1)] - \delta - x, \end{align}

where $\theta(x)$ is the Heaviside function. We will now step through a single oscillation to see how long it takes (the period) and how high $x$ gets (the amplitude). We will start at $\mathrm{d}x/\mathrm{d}t < 0$ and $x = 1$, i.e., on our way down in $x$. We will call this $t = 0$. No X will be produced until $t = \tau$, since until that time, $x(t-\tau) > 1$. This portion of the oscillation takes time $\tau$. During this time, the concentration of $x$ continues to drop, roughly to zero. So, now, we are at time $t = \tau$ with $x \approx 0$.

The next part of the oscillation involves production of X.

  • i) Show that the amount of time it takes to get to $x = 1$ is $t_1 = \ln[(\beta - \delta) / (\beta - \delta - 1)]$. Hint: Write down and solve the ODE for $x(t)$ in this regime in which $x$ is being produced at rate $\beta$ and degraded at rate $\delta + x$. Recall that the solution to the ODE $\dot{x} = a + bx$ is $x(t) = \left(a/b+x(0)\right) \mathrm{e}^{bt} - a/b$.
  • ii) From time $t = \tau + t_1$, X will continue to be produced until $t = 2\tau + t_1$, when production of X halts. Explain why this is the case.
  • iii) Show that $x(2\tau + t_1) = (\beta - \delta) + (1 - \beta + \delta) \mathrm{e}^{-\tau}$. This is the maximal value of $x$, and is therefore the amplitude. So, we define $A = (\beta - \delta) + (1 - \beta + \delta) \mathrm{e}^{-\tau}$.
  • iv) Coming down from the peak, with X no longer being produced, we complete the oscillation by returning to $x = 1$. Show that the time it takes to return to $x = 1$ from the peak of the oscillation is $t_2 = \ln[(A+\delta)/(1+\delta)]$.

  • v) The total time of the oscillation (the period) is therefore $T = 2\tau + t_1 + t_2$. Show that in the limit of $\beta,\delta \gg 1$,

    \begin{align} T \approx 2\tau + \ln\left[\frac{\beta}{\delta}\left(1 - \mathrm{e}^{-\tau}\right) + \mathrm{e}^{-\tau}\right]. \end{align}

  • vi) How do your analytical amplitude and period compare to what you computed in your numerical analysis?

e) These analytical results reveal a strategy for tuning the oscillator.

  • i) How could you tune the period and amplitude of the oscillator independently? Consider cases where $\tau \gg 1$ and $\tau \ll 1$.
  • ii) Change parameters to make the amplitude of the oscillator twice as large while keeping the same period as in part (c) for both parameter sets considered there. Numerically solve, plot the result, and explain your strategy.
  • iii) Change parameters to make the period of the oscillator twice as long while keeping the same amplitude as in part (c) for both parameter sets considered there. Numerically solve, plot the result, and explain your strategy.
  • iv) Based on your analysis, why is it important to have enzymatic degradation of X to have tunable oscillations?


5.2 Noise in a switchable promoter, 50 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$.

  • i) 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}
  • ii) 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.

  • i) 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.
  • ii) 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. 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}$.

  • i) 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.

  • ii) 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.
  • iii) Can you get qualitatively different steady state distributions of RNA copy number by varying $k_\mathrm{on}$ and $k_\mathrm{off}$? Plot representative distributions.