Homework 7 - Stochastic Operon Simulation - First Draft Due 1pm Friday October 12. Final Draft Due 5pm Wednesday, October 17

And Now with Extra Credit

We will follow McAdams and Arkin in their application of Gillespie's method to a simple model of gene expression. In particular we capture the transcription, translation and dimerization at gene 1 in figure 1 via the following reactions 

  R + PR -> RPo          % RNAP promoter binding and opening
  R + PR <- RPo          % RNAP promoter closing and unbinding
  RPo -> 10P + R + PR    % translation of protein and freeing of R and PR
  P + P -> D             % dimerization
  P + P <- D             % dedimerization
  P ->                   % protein degradation

We collect these variables into the vector 

   x = [R PR RPo P D]

and follow their evolution from the initial levels 

   x0 = [10 1 0 0 0]

and propensities 

   c = [2 1 4 2 0.5 0.05]

Working from the outside in, you will write a function called gillrunner that takes as arguments 

        tfin,   the duration of the simulation
        tinc,	the increment between interpolated time points
        nr,     the number of Gillespie runs
        rtab,	a cell array encoded the 6 reactions above
        x0,     the 5x1 vector above
        c,      the 6x1 vector above (or your variant)

and call your subfunction mygill nr times and plots (as below) the the mean and plus/minus one standard deviation of the dimer solutions.

Your mygill subfunction should behave like 

    function [t,dimer] = mygill(tfin, rtab, x, k)

and should call on its associated update function and should not require any if clauses.

In the absense of inhibition the dimer count is roughly increaing. To complete this assignment I ask you to augment your code to incorporate repression of PR by D via the two additional reactions 

  PR + D -> I            % D binds PR and so R can not read
  I -> PR + D            % D unbinds PR

I would start from I=0 and assume propensities of 0.01 and 0.005 and produce something like that below. (I believe I used nr=4. You are encouraged to experiment with propensities and run numbers.)

Your work will be graded as follows: 

        First draft
                2 pts for gilldriver loop
                4 pts for mygill outline using rtab cell and update

        Final draft 
                8 pts for header CONTAINING detailed USAGE
                8 pts for further comments in code 
                4 pts for indentation
                10 pts for correct gillrunner
                10 pts for correct computation of reaction number
                8 pts for correct coding of update

	        6 pts for two jpeg plots corresponding to the free dimer
                      production, and to its self-regulated production.

        EXTRA CREDIT - PLEDGED
                15 points for a 500 word essay. The essay must be typed and
                have a title, introduction and body. The backbone of the
                essay should be an overview of the paper by McAdams and Arkin. 
                Into that foundation you should weave detailed reference to 
                your own Boolean and Gillespien gene work and discuss its
                connection to the gene knock-out work that this week received
                the Nobel prize in Physiology or Medicine.
                You may neither give nor receive aid from other students. You should consult
                and cite both paper and electronic media. Your citation list
                does not count toward the 500 word minimum.