Code

From 2007.igem.org

(Difference between revisions)
 
Line 1: Line 1:
-
function f = lambda_circuit(tfin,tinc,nr,x0)
+
function f = lambda_circuit(tfin,tinc,nr,x0,off)
 +
% x0 = [1 0 0 0 0 1 0 0 0 0 0 0 0];
%
%
-
% f = lambda_circuit(tfin,tinc,nr,x0)
 
-
% Calculates the Gillespie realizations of the
 
-
% Rice IGEM phage construct.
 
-
%
 
-
% tfin - final time of simulation (min)
 
-
% tinc - sampling rate (min)
 
-
% nr - number of Gillespie runs
 
-
% initial conditions
 
% x = [P_tetR_Open P_tetR_Closed mRNA_tetR TetR TetR2 P_cI_open ...
% x = [P_tetR_Open P_tetR_Closed mRNA_tetR TetR TetR2 P_cI_open ...
%        P_cI_closed mRNA_cI CI]
%        P_cI_closed mRNA_cI CI]
-
% for example
 
%
%
-
% lambda_circuit(500,10,20,[1 0 0 0 0 1 0 0 0 0])
 
-
 
-
 
% x(1) = P_tetR_Open
% x(1) = P_tetR_Open
% x(2) = P_tetR_Closed
% x(2) = P_tetR_Closed
Line 26: Line 15:
% x(9) = CI
% x(9) = CI
% x(10) =  CI2
% x(10) =  CI2
 +
% x(11) = P_tet_closed
 +
% x(12) = Tc
 +
% x(13) = TetR2-Tc-TetR2
 +
% x(14) = Tc-TetR2_out
 +
% x(15) = Tc_out
 +
% x(16) = TetR2_out
 +
 +
 +
% Reactions
% Reactions
%
%
% 1) P_tetR_open + RNAP --> mRNA_tetR + RNAP
% 1) P_tetR_open + RNAP --> mRNA_tetR + RNAP
-
%  x(1) ---> x(1) + x(3)
+
%  x(1) ---> x(1) + x(3), 2/min
%
%
% 2) mRNA_tetR + RIB ---> TetR + RIB + mRNA_tetR
% 2) mRNA_tetR + RIB ---> TetR + RIB + mRNA_tetR
-
%  x(3) ---> x(3) + x(4)
+
%  x(3) ---> x(3) + x(4), 2/min
%
%
% 3) mRNA_tetR -->
% 3) mRNA_tetR -->
-
%  x(3) --->
+
%  x(3) --->, 1/min
%
%
% 4) 2TetR ---> TetR2
% 4) 2TetR ---> TetR2
-
%  2x(4) ---> x(5)
+
%  2x(4) ---> x(5), 0.03/nM/min
%
%
% 5) TetR2 ---> 2TetR
% 5) TetR2 ---> 2TetR
-
%  x(5) ---> 2x(4)
+
%  x(5) ---> 2x(4) 0.001/min
%
%
% 6) TetR --->
% 6) TetR --->
-
%  x(4) --->
+
%  x(4) ---> 0.1/min
%
%
% 7) TetR2 --->
% 7) TetR2 --->
-
%  x(5) ---->
+
%  x(5) ----> 0.1/min
%
%
% 8) TetR2 + p_cI_open ---> P_cI_closed
% 8) TetR2 + p_cI_open ---> P_cI_closed
-
%  x(5) + x(6) ---> x(7)
+
%  x(5) + x(6) ---> x(7), 0.2/nM/min
%
%
% 9) p_cI_closed ---> p_cI_open + TetR2
% 9) p_cI_closed ---> p_cI_open + TetR2
-
%  x(7) ---> x(5) + x(6)
+
%  x(7) ---> x(5) + x(6), 0.1/min
%
%
% 10) p_cI_closed + RNAP ---> p_cI_closed + RNAP + mRNA_cI
% 10) p_cI_closed + RNAP ---> p_cI_closed + RNAP + mRNA_cI
-
%  x(7) ---> x(7) + x(8)
+
%  x(7) ---> x(7) + x(8), 0.02/min
%
%
% 11) p_cI_open + RNAP ---> p_cI_open + RNAP + mRNA_cI
% 11) p_cI_open + RNAP ---> p_cI_open + RNAP + mRNA_cI
-
%  x(6) ---> x(6) + x(8)
+
%  x(6) ---> x(6) + x(8), 2/min
%
%
% 12) mRNA_cI --->
% 12) mRNA_cI --->
-
%  x(8) --->
+
%  x(8) --->, 1/min
%
%
% 13) mRNA_cI + RIB ---> mRNA_cI + CI + RIB
% 13) mRNA_cI + RIB ---> mRNA_cI + CI + RIB
-
%  x(8) ---> x(8) + x(9)
+
%  x(8) ---> x(8) + x(9), 2/min
%
%
% 14) 2CI ---> CI2
% 14) 2CI ---> CI2
-
%  2x(9) ---> x(10)
+
%  2x(9) ---> x(10), 0.03/nM/min
%
%
% 15) CI2 ---> 2CI
% 15) CI2 ---> 2CI
-
%  x(10) ---> 2x(9)
+
%  x(10) ---> 2x(9), 0.001/min
%
%
% 16) CI --->
% 16) CI --->
-
%  x(9) --->
+
%  x(9) ---> 0.1/min
%
%
% 17) CI2 --->
% 17) CI2 --->
-
%  x(10) --->
+
%  x(10) ---> 0.1/min
%
%
 +
% 18) P_tet_open + TetR2 ---> P_tet_closed
 +
%  x(1) + x(5) ---> x(11), 0.2/na M/min
 +
%
 +
% 19) P_tet_closed ---> P_tet_open + TetR2
 +
%  x(11) ---> x(1) + x(5), 0.1/min
 +
%
 +
% 20) Tc + 2 TetR2 ---> Tc-TetR2
 +
%  x(12) + 2x(5) ---> x(13), 0.05/nM/min
 +
%
 +
% 21) Tc-TetR2 ---> Tc + 2 TetR2
 +
%  x(13) ---> x(12) + 2 x(5), 0.1/min
 +
%
 +
% 22) Tc-TetR2 ---> Tc-TetR2_out
 +
%  x(13) ----> x(14)  0.1/min
 +
%
 +
% 23) P_tet_closed ---> mRNA_tetR + P_tet_closed
 +
%  x(11) ---> x(11) + x(3), 0.02/min
 +
%
 +
% 24) Tc_out ---> Tc
 +
%  x(15) ---> x(12), .01/min
 +
%
 +
% 25) Tc-TetR2_out --> Tc_out + 2 TetR2_out
 +
%  x(14) ---> x(15) + 2 x(16)
 +
%
 +
% 26) Tc_out + 2TetR2_out ---> Tc-TetR2_out
 +
%  X(15) + 2 x(16) ---> x(14)
 +
%
 +
% 26) Tc-TetR2 --->
 +
% 
 +
 +
close all;
-
% Reaction Rates
+
transo= 2;      %transcription from open promoter
-
k = [2 2 .1 0.003 .1 .005 .01 .2 .1 .002 2 2 .1 .003 .1 .01 .01];
+
tranl = 2;     %translation
-
%k = [2 2 .5 0.03 .1 .05 .05 .2 .1 .02 2 2 .5 .03 .1 .05 .05];
+
mdeg = .5;      % mRNA degredation
 +
dime = .03;    % dimerization
 +
udim = 0.001;  % undimerization
 +
pdeg = .05;    % protein degredation
 +
proc = 0.2;    % closed promoter formation rate
 +
proo = 0.1;    % open promoter formation rate
 +
transc = 0.02;  % transcription from closed promoter
 +
induc = 0.5;  % inducer binding
 +
uninduc = 0.1; % inducer unbinding
 +
k = [transo    % 1
 +
    tranl      % 2
 +
    mdeg        % 3
 +
    dime        % 4
 +
    udim        % 5
 +
    pdeg        % 6
 +
    pdeg        % 7
 +
    proc        % 8
 +
    proo        % 9
 +
    transc      % 10
 +
    transo      % 11
 +
    mdeg        % 12
 +
    tranl      % 13
 +
    dime        % 14
 +
    udim        % 15
 +
    pdeg        % 16
 +
    pdeg        % 17
 +
    proc/5    % 18 weaken TetR2 binding/unbinding to tetR promoter
 +
    proo/5    % 19
 +
    induc      % 20
 +
    uninduc    % 21
 +
    1        % 22 Extrusion of Tc-TetR out of cell
 +
    transc      % 23
 +
    .5          % 24 diffusion of Tc into cell
 +
    uninduc    % 25
 +
    induc]';    % 26
-
% stochiometric action of each reactions
+
k([14 15 17]) = [0 0 0];   
-
% rtab{i}(1:2:end) = species in reaction i
+
-
% rtab{i}(2:2:end} = action of species in reaction i
+
rtab = {[1 0 3 1]      % 1
rtab = {[1 0 3 1]      % 1
Line 106: Line 167:
     [9 2 10 -1]        % 15
     [9 2 10 -1]        % 15
     [9 -1]              % 16
     [9 -1]              % 16
-
     [10 -1]};           % 17
+
     [10 -1]            % 17
 +
    [1 -1 5 -1 11 1]    % 18
 +
    [1 1 5 1 11 -1]    % 19
 +
    [5 -2 12 -1 13 1]  % 20
 +
    [5 2 12 1 13 -1]    % 21
 +
    [13 -1 14 1]        % 22
 +
    [3 1 11 0]          % 23
 +
    [12 1 15 -1]        % 24
 +
    [14 -1 15 1 16 2]  % 25
 +
    [14 1 15 -1 16 -2]};% 26
-
% initialize sampling rate
 
tvec = 1:tinc:tfin;
tvec = 1:tinc:tfin;
Line 116: Line 185:
P0 = cell(1,L);
P0 = cell(1,L);
 +
global x_num
-
%initialize probabilities for x at sampled times
+
x_num = length(x0);
-
for i = 1:L
+
 
-
   
+
if off == 1,
-
    P0{i} = zeros(10,5000);
+
 
 +
    for i = 1:L
 +
 
 +
        P0{i} = zeros(x_num,5000);
 +
 
 +
    end
end
end
Line 130: Line 205:
     disp(['Proceeding with run ' num2str(j)])  
     disp(['Proceeding with run ' num2str(j)])  
-
     [t,xout] = gomcreg(tfin, rtab, x0, k);     %Gillespie step 
+
     [t,xout] = gomcreg(tfin, rtab, x0, k);      
      
      
-
     xinterp = interp0(t,xout,tvec);             %0-order interpolation
+
     xinterp = interp0(t,xout,tvec);
-
   
+
-
    X([1:10]+(j-1)*10,:) = xinterp;            %store sampled runs
+
      
      
-
     P0 = updateP(P0,xinterp,L);                 %update probabilites
+
     X([1:x_num]+(j-1)*x_num,:) = xinterp;
 +
   
 +
    if off == 1,
 +
        P0 = updateP(P0,xinterp,L);
 +
    end
 +
       
end
end
-
  m = zeros(10,length(tvec));
+
  m = zeros(x_num,length(tvec));
  s = m;
  s = m;
-
for j=1:10,
+
for j=1:x_num,
-
     m(j,:) = mean(X(j:10:end,:),1);  % compute mean
+
     m(j,:) = mean(X(j:x_num:end,:),1);  % compute ensemble mean
-
     s(j,:) = std(X(j:10:end,:),1);  % compute standard deviation
+
     s(j,:) = std(X(j:x_num:end,:),1);  % compute ensemble mean
end
end
-
slab = {'P_tetR_Open' 'P_tetR_Closed' 'mRNA_tetR' 'TetR' 'TetR2' ...
 
-
    'P_cI_open' 'P_cI_closed' 'mRNA_cI' 'CI' 'CI2'};
 
-
i=0;
 
 +
slab = {'P_{tetR}_{Open}' 'P_{tetR}_{Closed}' 'mRNA_{tetR}' 'TetR' 'TetR_2' ...
 +
    'P_{cI}_{open}' 'P_{cI}_{closed}' 'mRNA_{cI}' 'CI' 'CI_2' 'P_{tetR}_{closed}' 'Tc' 'TetR_2-Tc' 'Tc-TetR_{out}' 'Tc_{out}' 'TetR_{out}'};
-
%Graph mean and std for the runs
+
i=0;
-
for j = [3 4 5 8 9 10]
+
for j = [ 4 5 9 12 13 14 15 16]
     i=i+1;
     i=i+1;
-
     subplot(2,3,i)
+
     subplot(2,4,i)
     plot(tvec,m(j,:),'b')
     plot(tvec,m(j,:),'b')
     hold on
     hold on
Line 171: Line 248:
end
end
-
pause
+
if off == 1,
 +
 
 +
 
 +
    for t = 1:L
 +
 
 +
        A = P0{t};
 +
 
 +
        p0 = (sum(A'));
 +
 
 +
        p0 = 1./p0;
 +
 
 +
        for i = 1:x_num
 +
 
 +
            P0{t}(i,:) = P0{t}(i,:)*p0(i);
 +
 
 +
        end
 +
 
 +
    end
 +
 
 +
    i =1;
 +
    for t = 1:10:L
 +
 
 +
        i = i+1;
 +
 
 +
        figure(i)
 +
 
 +
        graphP(P0{t},m(:,t)+5*s(:,t),tvec(t));
-
%recover and plot the probabilites
 
-
for t = 1:L
 
-
   
 
-
    A = P0{t};
 
-
 
 
-
    p0 = (sum(A'));
 
-
   
 
-
    p0 = 1./p0;
 
-
   
 
-
    for i = 1:10
 
-
       
 
-
        P0{t}(i,:) = P0{t}(i,:)*p0(i);
 
-
       
 
     end
     end
-
   
 
-
end
 
-
i =1;
 
-
for t = 1:10:L
 
-
   
 
-
    i = i+1;
 
-
   
 
-
    figure(i)
 
-
   
 
-
    graphP(P0{t},m(:,t)+5*s(:,t));
 
-
   
 
-
   
 
-
   
 
end
end
 +
f = m(:,end);
 +
 +
%keyboard
return
return
Line 219: Line 299:
function [t,xout] = gomcreg(tfin, rtab, x, c)
function [t,xout] = gomcreg(tfin, rtab, x, c)
-
t = zeros(1,1e6);  % preallocate space for t and dimer, for speed
+
global x_num
-
xout = zeros(10,1e6);  
+
 
 +
t = zeros(1,1e7);  % preallocate space for t and dimer, for speed
 +
xout = zeros(x_num,1e7);  
t(1) = 0;
t(1) = 0;
xout(:,1) = x';
xout(:,1) = x';
Line 231: Line 313:
     a = c.*h; % reaction probabilities
     a = c.*h; % reaction probabilities
     a0 = sum(a);
     a0 = sum(a);
 +
     
 +
%    if mod(i,10000) == 0,
 +
%        disp(['a0= ',num2str(a0), 'on run ', num2str(i)]);
 +
%    end
 +
   
      
      
-
    %if mod(i,1000) == 0,
 
-
    %    disp(['a0= ',num2str(a0), 'on run ', num2str(i)]);
 
-
    %end
 
      
      
     r1 = rand;
     r1 = rand;
Line 250: Line 334:
      
      
     x(metabs) = x(metabs) + action;
     x(metabs) = x(metabs) + action;
-
 
+
   
-
     xout(:,i) = x';
+
     xout(:,i) = x'; % record dimer count
end
end
Line 257: Line 341:
disp(['    took ' num2str(i) ' Gillespie steps'])  % talk to user
disp(['    took ' num2str(i) ' Gillespie steps'])  % talk to user
-
t = t(1:i);
+
t = t(1:i); % only return the computed t and dimer vals
xout = xout(:,1:i);
xout = xout(:,1:i);
Line 277: Line 361:
h(9) = x(7);
h(9) = x(7);
h(10) = x(7);
h(10) = x(7);
-
h(11) = x(7);
+
h(11) = x(6);
h(12) = x(8);
h(12) = x(8);
h(13) = x(8);
h(13) = x(8);
Line 284: Line 368:
h(16) = x(9);
h(16) = x(9);
h(17) = x(10);
h(17) = x(10);
 +
h(18) = x(1)*x(5);
 +
h(19) = x(11);
 +
h(20) = x(12)*x(5)*(x(5)-1)/2;
 +
h(21) = x(13);
 +
h(22) = x(13);
 +
h(23) = x(11);
 +
h(24) = x(15);
 +
h(25) = x(14);
 +
h(26) = x(15)*(x(16)-1)/2;
function P = updateP(P0,xout,L)
function P = updateP(P0,xout,L)
 +
 +
global x_num
for t = 1:L
for t = 1:L
      
      
-
     for i = 1:10
+
     for i = 1:x_num
       P0{t}(i,xout(i,t)+1) = P0{t}(i,xout(i,t)+1) + 1;
       P0{t}(i,xout(i,t)+1) = P0{t}(i,xout(i,t)+1) + 1;
Line 333: Line 428:
-
function f = graphP(P,v)
+
function graphP(P,v,t)
[m n] = size(P);
[m n] = size(P);
-
slab = {'P_tetR_Open' 'P_tetR_Closed' 'mRNA_tetR' 'TetR' 'TetR2' ...
+
slab = {'P_{tetR}_{Open}' 'P_{tetR}_{Closed}' 'mRNA_{tetR}' 'TetR' 'TetR_2' ...
-
     'P_cI_open' 'P_cI_closed' 'mRNA_cI' 'CI' 'CI2'};
+
     'P_{cI}_{open}' 'P_{cI}_{closed}' 'mRNA_{cI}' 'CI' 'CI_2' 'P_{tetR}_{closed}' 'Tc' 'TetR_2-Tc' 'Tc-TetR_out' 'Tc_out' 'TetR_out'};
 +
 
i=0;
i=0;
-
for j = [3 4 5 8 9 10]
+
for j = [ 4 5 9 10 12 13]
     i=i+1;
     i=i+1;
     subplot(2,3,i)
     subplot(2,3,i)
Line 352: Line 448:
     ax(2) = v(j)+1;
     ax(2) = v(j)+1;
     axis(ax)
     axis(ax)
 +
    title([num2str(t),' (min)']);
     hold off
     hold off
end
end

Latest revision as of 03:46, 27 October 2007

function f = lambda_circuit(tfin,tinc,nr,x0,off) % x0 = [1 0 0 0 0 1 0 0 0 0 0 0 0]; % % x = [P_tetR_Open P_tetR_Closed mRNA_tetR TetR TetR2 P_cI_open ... % P_cI_closed mRNA_cI CI] % % x(1) = P_tetR_Open % x(2) = P_tetR_Closed % x(3) = mRNA _tetR % x(4) = TetR % x(5) = TetR2 % x(6) = P_cI_open % x(7) = p_cI_closed % x(8) = mRNA_cI % x(9) = CI % x(10) = CI2 % x(11) = P_tet_closed % x(12) = Tc % x(13) = TetR2-Tc-TetR2 % x(14) = Tc-TetR2_out % x(15) = Tc_out % x(16) = TetR2_out



% Reactions % % 1) P_tetR_open + RNAP --> mRNA_tetR + RNAP % x(1) ---> x(1) + x(3), 2/min % % 2) mRNA_tetR + RIB ---> TetR + RIB + mRNA_tetR % x(3) ---> x(3) + x(4), 2/min % % 3) mRNA_tetR --> % x(3) --->, 1/min % % 4) 2TetR ---> TetR2 % 2x(4) ---> x(5), 0.03/nM/min % % 5) TetR2 ---> 2TetR % x(5) ---> 2x(4) 0.001/min % % 6) TetR ---> % x(4) ---> 0.1/min % % 7) TetR2 ---> % x(5) ----> 0.1/min % % 8) TetR2 + p_cI_open ---> P_cI_closed % x(5) + x(6) ---> x(7), 0.2/nM/min % % 9) p_cI_closed ---> p_cI_open + TetR2 % x(7) ---> x(5) + x(6), 0.1/min % % 10) p_cI_closed + RNAP ---> p_cI_closed + RNAP + mRNA_cI % x(7) ---> x(7) + x(8), 0.02/min % % 11) p_cI_open + RNAP ---> p_cI_open + RNAP + mRNA_cI % x(6) ---> x(6) + x(8), 2/min % % 12) mRNA_cI ---> % x(8) --->, 1/min % % 13) mRNA_cI + RIB ---> mRNA_cI + CI + RIB % x(8) ---> x(8) + x(9), 2/min % % 14) 2CI ---> CI2 % 2x(9) ---> x(10), 0.03/nM/min % % 15) CI2 ---> 2CI % x(10) ---> 2x(9), 0.001/min % % 16) CI ---> % x(9) ---> 0.1/min % % 17) CI2 ---> % x(10) ---> 0.1/min % % 18) P_tet_open + TetR2 ---> P_tet_closed % x(1) + x(5) ---> x(11), 0.2/na M/min % % 19) P_tet_closed ---> P_tet_open + TetR2 % x(11) ---> x(1) + x(5), 0.1/min % % 20) Tc + 2 TetR2 ---> Tc-TetR2 % x(12) + 2x(5) ---> x(13), 0.05/nM/min % % 21) Tc-TetR2 ---> Tc + 2 TetR2 % x(13) ---> x(12) + 2 x(5), 0.1/min % % 22) Tc-TetR2 ---> Tc-TetR2_out % x(13) ----> x(14) 0.1/min % % 23) P_tet_closed ---> mRNA_tetR + P_tet_closed % x(11) ---> x(11) + x(3), 0.02/min % % 24) Tc_out ---> Tc % x(15) ---> x(12), .01/min % % 25) Tc-TetR2_out --> Tc_out + 2 TetR2_out % x(14) ---> x(15) + 2 x(16) % % 26) Tc_out + 2TetR2_out ---> Tc-TetR2_out % X(15) + 2 x(16) ---> x(14) % % 26) Tc-TetR2 ---> %

close all;

transo= 2;  %transcription from open promoter tranl = 2;  %translation mdeg = .5;  % mRNA degredation dime = .03;  % dimerization udim = 0.001;  % undimerization pdeg = .05;  % protein degredation proc = 0.2;  % closed promoter formation rate proo = 0.1;  % open promoter formation rate transc = 0.02;  % transcription from closed promoter induc = 0.5;  % inducer binding uninduc = 0.1;  % inducer unbinding

k = [transo  % 1 

   tranl       % 2
   mdeg        % 3
   dime        % 4
   udim        % 5
   pdeg        % 6
   pdeg        % 7
   proc        % 8
   proo        % 9
   transc      % 10
   transo      % 11
   mdeg        % 12
   tranl       % 13
   dime        % 14
   udim        % 15
   pdeg        % 16
   pdeg        % 17
   proc/5     % 18 weaken TetR2 binding/unbinding to tetR promoter
   proo/5     % 19
   induc       % 20
   uninduc     % 21
   1        % 22 Extrusion of Tc-TetR out of cell
   transc      % 23
   .5          % 24 diffusion of Tc into cell
   uninduc     % 25
   induc]';    % 26

k([14 15 17]) = [0 0 0];

rtab = {[1 0 3 1]  % 1 

   [3 0 4 1]           % 2
   [3 -1]              % 3
   [4 -2 5 1]          % 4
   [4 2 5 -1]          % 5  
   [4 -1]              % 6
   [5 -1]              % 7
   [5 -1 6 -1 7 1]     % 8
   [5 1 6 1 7 -1]      % 9
   [7 0 8 1]           % 10
   [6 0 8 1]           % 11
   [8 -1]              % 12
   [8 0 9 1]           % 13
   [9 -2 10 1]         % 14
   [9 2 10 -1]         % 15
   [9 -1]              % 16
   [10 -1]             % 17
   [1 -1 5 -1 11 1]    % 18
   [1 1 5 1 11 -1]     % 19
   [5 -2 12 -1 13 1]   % 20
   [5 2 12 1 13 -1]    % 21
   [13 -1 14 1]        % 22
   [3 1 11 0]          % 23
   [12 1 15 -1]        % 24
   [14 -1 15 1 16 2]   % 25
   [14 1 15 -1 16 -2]};% 26


tvec = 1:tinc:tfin;

L = ceil(tfin/tinc);

P0 = cell(1,L);

global x_num

x_num = length(x0);

if off == 1, 

   for i = 1:L

       P0{i} = zeros(x_num,5000);

   end

end


for j=1:nr 

   disp(['Proceeding with run ' num2str(j)]) 

   [t,xout] = gomcreg(tfin, rtab, x0, k);        
   
   xinterp = interp0(t,xout,tvec);
   
   X([1:x_num]+(j-1)*x_num,:) = xinterp;
   
   if off == 1,
       P0 = updateP(P0,xinterp,L);
   end


end 

m = zeros(x_num,length(tvec));
s = m;

for j=1:x_num, 

    m(j,:) = mean(X(j:x_num:end,:),1);  % compute ensemble mean
    s(j,:) = std(X(j:x_num:end,:),1);  % compute ensemble mean

end


slab = {'P_{tetR}_{Open}' 'P_{tetR}_{Closed}' 'mRNA_{tetR}' 'TetR' 'TetR_2' ... 

   'P_{cI}_{open}' 'P_{cI}_{closed}' 'mRNA_{cI}' 'CI' 'CI_2' 'P_{tetR}_{closed}' 'Tc' 'TetR_2-Tc' 'Tc-TetR_{out}' 'Tc_{out}' 'TetR_{out}'};

i=0;

for j = [ 4 5 9 12 13 14 15 16] 

   i=i+1;
   subplot(2,4,i)
   plot(tvec,m(j,:),'b')
   hold on
   plot(tvec,m(j,:)+s(j,:),'r')
   plot(tvec,m(j,:)-s(j,:),'r')
   xlabel('t')
   ylabel(slab{j})
   ax = axis;
   ax(2) = tfin;
   axis(ax)
   hold off

end

if off == 1, 


   for t = 1:L

       A = P0{t};

       p0 = (sum(A'));

       p0 = 1./p0;

       for i = 1:x_num

           P0{t}(i,:) = P0{t}(i,:)*p0(i);

       end

   end

   i =1;
   for t = 1:10:L

       i = i+1;

       figure(i)

       graphP(P0{t},m(:,t)+5*s(:,t),tvec(t));



   end

end

f = m(:,end);

%keyboard

return

% % straight Gillespie on the self-regulating gene % % usage [t,dimer] = gomcreg(tfin, rtab, x, c) % % where tfin, rtab, x, and k are as described above % and t and dimer are the resulting time and dimer count %

function [t,xout] = gomcreg(tfin, rtab, x, c)

global x_num

t = zeros(1,1e7);  % preallocate space for t and dimer, for speed xout = zeros(x_num,1e7); t(1) = 0; xout(:,1) = x';

i = 1;

while t(i) < tfin, % proceed until tfin 

   h = update(x);              % # of distint reactant combinations
   a = c.*h;			% reaction probabilities
   a0 = sum(a);

% if mod(i,10000) == 0, % disp(['a0= ',num2str(a0), 'on run ', num2str(i)]); % end 


   r1 = rand;
   i = i + 1;
   tau = log(1/r1)/a0;		% the time to next reaction
   t(i) = t(i-1) + tau;	% record time

   r2 = rand;
   b = cumsum(a)/a0;
   reac = find(b>r2,1);		% the reaction number

   metabs = rtab{reac}(1:2:end);       % the associated reactants
   action = rtab{reac}(2:2:end);       % the associated action
   
   
   x(metabs) = x(metabs) + action;
   
   xout(:,i) = x';		% record dimer count

end

disp([' took ' num2str(i) ' Gillespie steps'])  % talk to user

t = t(1:i); % only return the computed t and dimer vals xout = xout(:,1:i);

return



function h = update(x) h(1) = x(1); h(2) = x(3); h(3) = x(3); h(4) = x(4)*(x(4)-1)/2; h(5) = x(5); h(6) = x(4); h(7) = x(5); h(8) = x(5)*x(6); h(9) = x(7); h(10) = x(7); h(11) = x(6); h(12) = x(8); h(13) = x(8); h(14) = x(9)*(x(9)-1)/2; h(15) = x(10); h(16) = x(9); h(17) = x(10); h(18) = x(1)*x(5); h(19) = x(11); h(20) = x(12)*x(5)*(x(5)-1)/2; h(21) = x(13); h(22) = x(13); h(23) = x(11); h(24) = x(15); h(25) = x(14); h(26) = x(15)*(x(16)-1)/2;


function P = updateP(P0,xout,L)


global x_num

for t = 1:L 

   for i = 1:x_num

      P0{t}(i,xout(i,t)+1) = P0{t}(i,xout(i,t)+1) + 1;
      
   end

end

P = P0;

return

function yi = interp0(x,y,xi)

[m n] = size(y);

yi = zeros(m,length(xi));

for i = 1:m, 

   for j = 1:length(xi)
       
       index = find( xi(j) <= x);

       if index(1) == 1,
       
           yi(i,j) = y(i,index(1));
       
       else

           yi(i,j) = y(i,index(1)-1);
   
       end

   end

end

return


function graphP(P,v,t)

[m n] = size(P);

slab = {'P_{tetR}_{Open}' 'P_{tetR}_{Closed}' 'mRNA_{tetR}' 'TetR' 'TetR_2' ... 

   'P_{cI}_{open}' 'P_{cI}_{closed}' 'mRNA_{cI}' 'CI' 'CI_2' 'P_{tetR}_{closed}' 'Tc' 'TetR_2-Tc' 'Tc-TetR_out' 'Tc_out' 'TetR_out'};


i=0;

for j = [ 4 5 9 10 12 13] 

   i=i+1;
   subplot(2,3,i)
   plot(1:n,P(j,:),'b')
   hold on
   xlabel(slab(j))
   ylabel('P')
   ax = axis;
   ax(2) = v(j)+1;
   axis(ax)
   title([num2str(t),' (min)']);
   hold off

end

Retrieved from "http://2007.igem.org/wiki/index.php/Code"