In [1]:
# RUN ME FIRST!!!
from IPython.core.display import HTML
def css_styling():
    styles = open("custom.css", "r").read()
    return HTML(styles)
css_styling()
Out[1]:
/* From Lorena Barba's AeroPython course: https://github.com/barbagroup/AeroPython */

Homework 3

Part I Problem:

  1. Consider the following model for positive feedback gene regulation by a protein monomer: \begin{align*} &DNA_S + P \underset{\gamma_{S}}{\stackrel{\gamma_{F}}{\rightleftharpoons}} DNA_F & mRNA &\stackrel{\beta}{\to} mRNA + P\\ &DNA_S \stackrel{\alpha_S}{\to} DNA_S + mRNA & mRNA &\stackrel{\delta}{\to} \varnothing \\ &DNA_F \stackrel{\alpha_F}{\to} DNA_F + mRNA & P &\stackrel{\mu}{\to} \varnothing. \end{align*} Let $F(t)$ denote the probability the gene is in the fast transciption state, $DNA_F$; $S(t)$ is the probability the gene is in the slow transciption state, $DNA_S$; $m(t)$ is the amount of $mRNA$; and $P(t)$ is the amount of protein. Assume $\alpha_F > \alpha_S$.
    1. Write out the full set of ODEs for $\{S(t),F(t),m(t),P(t)\}$ corresponding to these reactions.
    2. Eliminate the equation for $F(t)$ using that $S(t) + F(t) = 1$.
    3. Apply a quasi-steady state approximation by setting $dS/dt$ to zero. What are the resulting equations for $dm/dt$ and $dP/dt$?
    4. Now set $dm/dt = 0$ too, what equation do you get for $dP/dt$?
    5. Show by graphical analysis as in class that there is only one non-negative steady state for $p(t)$. Make sure to justify why the curves have the shapes you draw.

Part II Problems:

  1. The toggle switch is a simple gene network where gene 1 can inhibit gene 2, and gene 2 can inhibit gene 1. Let $m_{1}(t)$ and $m_{2}(t)$ denote the number of mRNAs of species 1 and 2 respectively. Similarly, let $P_{1}(t)$ and $P_{2}(t)$ denote the numbers of proteins of each species. If we assume the promoter states are in quasi-steady state as in class, we get the following ODEs \begin{align} \frac{dm_{1}}{dt} &= \frac{ \alpha_{1}}{1 + (P_{2} / K_{1})^{n}} - \delta_{1} m_{1} + \alpha^{L}_{1},\\ \frac{dm_{2}}{dt} &= \frac{ \alpha_{2}}{1 + (P_{1} / K_{2})^{n}} - \delta_{2} m_{2} + \alpha^{L}_{2},\\ \frac{dP_{1}}{dt} &= \beta_{1} m_{1} - \mu_{1} P_{1},\\ \frac{dP_{2}}{dt} &= \beta_{2} m_{2} - \mu_{2} P_{2}. \end{align} The code below solves these ODEs using SciPy's `odeint` routine (the function defining the right hand side of the ODEs is `toggleF`). (`toggleProtOnlyF` returns the right hand side when $dm_1/dt=0$ and $dm_2/dt=0$ too. It is only used for plotting nullclines in $P_1$ vs $P_2$ space.)
    1. Write down a full set of reactions, including binding reactions for the repressor to the DNA promoter state, that give the preceding system of ODEs in the cases that 1) $n=1$ and 2) $n=2$. What, biologically, is the difference between using $n=1$ vs $n=2$? (You do not need to write the reaction rates, just the chemical reactions.)
    2. For $n=1$ solve the ODEs for the initial conditions in the code that `y0 = [0,0,0,0]`, `y0 = [0,0,20,0]`, and `y0 = [0,0,0,20]`. Now rerun each case for $n=2$. Use the six plots of $P_{1}(t)$ and $P_{2}(t)$ that are generated to explain the differences in behavior between the two cases. What happens at steady state? To help in explaining this, set `plotSS = True` in the code and run one of the $n=1$ cases and one of the $n=2$ cases. The new figure plots the two nullclines $dP_{1}/dt=0$ and $dP_{2}/dt=0$. Explain how they help determine the observed steady-state behavior. What do these results say about negative feedback by monomers vs. dimers?
    3. If we add a third gene, so that gene 1 represses gene 2, gene 2 represses gene 3, and gene 3 represses gene 1 we get the repressilator discussed in class (when $n=2$). Modify the code below to solve the repressilator system. (You may keep the same values for the parameters.) Show that if $P_{1}(0)=20$ and all the other variables are initially zero that you get a system that oscillates. Plot the three protein values vs. time on one figure, and turn in the figure. Why does the system oscillate? What happens if all the variables are initially zero? Why?
    4. How does changing the protein degradation rate in the repressilator influence the results? What about increasing the leakage rate to `5e-3`? Explain, biologically, what is happening in each of these cases, using figures from the simulations to justify your conclusions.

Note: I recommend leaving the code below as is and using the two code blocks below it for parts 2B, 2C, and 2D. Also, for 2C and 2D you might find it helpful to play around with different values of tf to understand what is going on in the model as you change parameters.

In [1]:
# RUN ME NEXT!!!
# we'll use the standalone plot interface to enable downloading/saving plots
# remember, plots may appear behind your web browser
%matplotlib 
Using matplotlib backend: MacOSX
In [2]:
# Toggle Switch Code for 2B, I'd recommend copying to a new code block before modifying.
import numpy as np
from matplotlib import pyplot as plt
import matplotlib as mpl
import scipy as sp
import scipy.integrate

def toggleF( y, t, alpha, alphaL, beta, delta, mu, K, n ):
    # RHS for two-gene toggle switch model.
    #   Parameters based on the repressilator.
    #   See Elowitz et al., Nature 2000.
    #   10/1/11 by Prof. Samuel Isaacson, isaacson@math.bu.edu
    
    # we assume y[0]=m1, y[1]=m2, y[2]=P1, y[3]=P2
    # f[0] = dm1/dt, f[1] = dm2/dt, f[2] = dP1/dt, f[3]=dP2/dt
    f    = np.zeros( np.shape(y) )
    f[0] = alpha[0] / (1 + (y[3] / K[0])**n[0]) - delta[0] * y[0] + alphaL[0]
    f[1] = alpha[1] / (1 + (y[2] / K[1])**n[1]) - delta[1] * y[1] + alphaL[1]
    f[2] = beta[0] * y[0] - mu[0] * y[2]
    f[3] = beta[1] * y[1] - mu[1] * y[3]
    
    return f


def toggleProtOnlyF( y, t, alpha, alphaL, beta, delta, mu, K, n ):
    # RHS for two-gene toggle switch model with mRNA dynamics at SS.
    #   Parameters based on the repressilator.
    #   See Elowitz et al., Nature 2000.
    # 10/2/15 by Prof. Samuel Isaacson, isaacson@math.bu.edu
        
    # mRNA steady states
    mss    = np.zeros( np.shape(y) )
    mss[0] = (alpha[0] / delta[0]) / (1 + (y[1] / K[0])**n[0]) + alphaL[0] / delta[0]
    mss[1] = (alpha[1] / delta[1]) / (1 + (y[0] / K[1])**n[1]) + alphaL[1] / delta[1]
    
    # f[0]=dP1/dt and f[1]=dP2/dt, given mRNA steady states
    f    = np.zeros( np.shape(y) )
    f[0] = beta[0] * mss[0] - mu[0] * y[0]
    f[1] = beta[1] * mss[1] - mu[1] * y[1]
    
    # note we could've written the previous three lines in a vectorized manner as:
    # f = beta * mss - mu * y
    
    return f


# PARAMETERS:

# the sum of the next two gives the transcription rate when there is no repression.
alpha  = np.array([.5,.5])                   # max repression transcription rate (per second)
alphaL = np.array([5e-4,5e-4])               # leak transcription rate (per second
delta  = np.log(2.) / np.array([120.,120.])  # mRNA degradation rates (per second)
beta   = 20. * delta                         # tranlation rates (per second), 20 proteins on average per mRNA life
mu     = np.log(2.) / np.array([600.,600.])  # protein degradation rate (per second)
K      = np.array([40.,40.])                 # dissociation constants for promoter regulation (molecule per sec)
n      = np.array([1.,1.])                   # Hill coefficient - 1 or 2
y0     = np.array([0.,0.,0.,20.])            # initial amounts y[0]=m1, y[1]=m2, y[2]=P1, y[3]=P2
tf     = 100000.                             # final time
plotSS = True                                # plot the nullclines dP1/dt=0 and dP2

############# END PARAMETERS


# times to get the solution at
t = np.linspace(0., tf, int(np.ceil(tf)) + 1)

# this defines a function f(y,t) corresponding to the derivative y'=f(y,t)
# we are essentially wrapping the function toggleF, which takes a bunch
# of extra parameters, to make a function that takes only y and t
f = lambda y,t : toggleF(y, t, alpha, alphaL, beta, delta, mu, K, n)

# Solve the ODEs using SciPy
y = sp.integrate.odeint(f, y0, t, atol=1e-12, rtol=1e-12)

# PLOTTING:

# pretty up plots
mpl.rcParams.update({'font.size': 12})  # this sets the default font size in plots
mpl.rcParams['axes.labelsize']  = 14    # make the xlabel/ylabel text size a bit bigger
mpl.rcParams['lines.linewidth'] = 2     # change the default linewidth in plots
    
# plotting
plt.figure()
plt.plot(t, y[:,0], 'k')
plt.plot(t, y[:,1], 'b')
plt.legend(('m1','m2'))
plt.xlabel('time (seconds)')
plt.ylabel('number of molecules')
plt.title('mRNA dynamics')

plt.figure()
plt.plot(t, y[:,2], 'k')
plt.plot(t, y[:,3], 'b')
plt.legend(('P1','P2'))
plt.xlabel('time (seconds)')
plt.ylabel('number of molecules')
plt.title('protein dynamics')

# nullcline plotting:
if plotSS:
    P1,P2 = np.meshgrid( 1.+np.linspace(0.,10000.,10000//25), 1.+np.linspace(0.,10000.,10000//25) )
    sz = np.shape( P1 )
    D1 = np.zeros( sz )
    D2 = np.zeros( sz )
    for i in range(0, sz[0]):
        for j in range(0, sz[1]):
            f = toggleProtOnlyF( np.array((P1[i,j], P2[i,j])), 1., alpha, alphaL, beta, delta, mu, K, n )
            D1[i,j] = f[0]
            D2[i,j] = f[1]

    fig = plt.figure()
    plt.xscale('log')
    plt.yscale('log')
    plt.contour(P1, P2, D1, [0.], colors='r')
    plt.contour(P1, P2, D2, [0.], colors='b')
    plt.xlabel('P1 steady state')
    plt.ylabel('P2 steady state')
    plt.title('P1 and P2 nullclines')    
In [4]:
# PUT YOUR CODE FOR 2B HERE
In [5]:
# PUT YOUR CODE FOR 2C HERE
  1. Expand the following function through all third derivative terms in a Taylor Series \begin{align*} f(t,y) &= \sin \left( t^{2} + y \right), &\text{about } t = 0, y = 0. \end{align*}
  2. Derive the third-order Taylor Series method for $f(t,y) = 2y-y^{2}$. Modify one of the Euler's method codes from the first homework or second lab to use the third-order Taylor series method. Then repeat part 3c) from the first homework.
  3. Show the following Runge-Kutta method is second order by Taylor Series expanding the associated truncation error about $t_{n}$: \begin{equation*} y_{n+1} = y_{n} + \Delta t \, f \left( t_{n}+\frac{\Delta t}{2}, y_{n} + \frac{\Delta t}{2} f(t_{n}, y_{n})\right). \end{equation*} (Hint, use two-dimensional Taylor series in expanding $f$.)
In [ ]: