Homework 8 - Deterministic vs. Stochastic Operon Simulation - First Draft Due 1pm Friday October 19. Final Draft Due 5pm Monday, October 22

We continue last week's investigation of a self-governing gene. In particular, we recall the following 8 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
  PR + D -> I            % D binds PR and so R can not read
  I -> PR + D            % D unbinds PR

involving the 6 species

   x = [R PR RPo P D I]

and propensities c(1) throught c(8).

Your first task is to write down the associated 6 ordinary differential equations (in the order they are listed in x above). I will give you the first one,

   d[R]/dt = -c(1)[R][PR] + c(2)[RPo] + c(3)[RPo]

Your next task is to augment last week's code so that it solves (via ode23) this ode system and compares the mean (stochastic) dimer level to the ode prediction. For example, using the initial counts and propensities

     x0 = [10 1 0 0 0 0]  and propensities   c = [2 1 4 2 0.5 0.05 .01 .005]

and nr=10, I arrive at the figure below

This result is promising, in that the dimer predictions appear to coincide. It would be nice to dig a little deeper and ask how the other reactants are fairing. At the stochastic level we know that each reactant may only take integer jumps and that R, PR, RPo and I in fact jump back and forth between 0 and 1. As taking ensemble means introduces intermediate values we instead compare our ode solution to single runss of Gillespie. The figure below was obtained with the same initial values and propemsities as above.

This figures cries out for explication. It shows early switching in R, PR and RPo, associated with stable early growth in D, despite the high volatility in P. And finally, once D achieves some critical level we see I switch on and stay on.

You will accomplish these two figures by modifying last week's gillrunner to behave like

    function mca2way(tfin, tinc, nr, x0, c)

and that produces a likeness of our first figure above when nr is greater than 1, and a likeness of our second figure (via subplot) when nr=1.

This function will call ode23 which will in turn call your subfunction

    function dx = mcaode(t,x,c)

Your work will be graded as follows:

        First draft
                3 pts for restructuring gillrunner
                3 pts for structure of mcaode

        Final draft 
                6 pts for header CONTAINING detailed USAGE
                  header of mcaode should give full system in detail
                6 pts for further comments in code
                4 pts for indentation
                12 pts for correct mca2way
                12 pts for correct mcaode

	        10 pts for two jpeg plots similar to those above

                4 pts for a typed, one paragraph contrast of the
                      stochastic and deterministic approaches to chemical
                      kinetics. Take it personally. Which method do
                      you prefer, and why?