From 2007.igem.org

function f = lambda_circuit(tfin,tinc,nr,x0) % % 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 ... % 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(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

% Reactions % % 1) P_tetR_open + RNAP --> mRNA_tetR + RNAP % x(1) ---> x(1) + x(3) % % 2) mRNA_tetR + RIB ---> TetR + RIB + mRNA_tetR % x(3) ---> x(3) + x(4) % % 3) mRNA_tetR --> % x(3) ---> % % 4) 2TetR ---> TetR2 % 2x(4) ---> x(5) % % 5) TetR2 ---> 2TetR % x(5) ---> 2x(4) % % 6) TetR ---> % x(4) ---> % % 7) TetR2 ---> % x(5) ----> % % 8) TetR2 + p_cI_open ---> P_cI_closed % x(5) + x(6) ---> x(7) % % 9) p_cI_closed ---> p_cI_open + TetR2 % x(7) ---> x(5) + x(6) % % 10) p_cI_closed + RNAP ---> p_cI_closed + RNAP + mRNA_cI % x(7) ---> x(7) + x(8) % % 11) p_cI_open + RNAP ---> p_cI_open + RNAP + mRNA_cI % x(6) ---> x(6) + x(8) % % 12) mRNA_cI ---> % x(8) ---> % % 13) mRNA_cI + RIB ---> mRNA_cI + CI + RIB % x(8) ---> x(8) + x(9) % % 14) 2CI ---> CI2 % 2x(9) ---> x(10) % % 15) CI2 ---> 2CI % x(10) ---> 2x(9) % % 16) CI ---> % x(9) ---> % % 17) CI2 ---> % x(10) ---> %

% Reaction Rates k = [2 2 .1 0.003 .1 .005 .01 .2 .1 .002 2 2 .1 .003 .1 .01 .01]; %k = [2 2 .5 0.03 .1 .05 .05 .2 .1 .02 2 2 .5 .03 .1 .05 .05];

% stochiometric action of each reactions % 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

   [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

% initialize sampling rate tvec = 1:tinc:tfin;

L = ceil(tfin/tinc);

P0 = cell(1,L);

%initialize probabilities for x at sampled times for i = 1:L

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

end

for j=1:nr

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

   [t,xout] = gomcreg(tfin, rtab, x0, k);      %Gillespie step  
   
   xinterp = interp0(t,xout,tvec);             %0-order interpolation
    
   X([1:10]+(j-1)*10,:) = xinterp;             %store sampled runs
   
   P0 = updateP(P0,xinterp,L);                 %update probabilites

end

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

for j=1:10,

    m(j,:) = mean(X(j:10:end,:),1);  % compute mean
    s(j,:) = std(X(j:10:end,:),1);  % compute standard deviation

end

slab = {'P_tetR_Open' 'P_tetR_Closed' 'mRNA_tetR' 'TetR' 'TetR2' ...

   'P_cI_open' 'P_cI_closed' 'mRNA_cI' 'CI' 'CI2'};

i=0;

%Graph mean and std for the runs

for j = [3 4 5 8 9 10]

   i=i+1;
   subplot(2,3,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

pause

%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

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

   i = i+1;
   
   figure(i)
   
   graphP(P0{t},m(:,t)+5*s(:,t));
   
   
   

end

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)

t = zeros(1,1e6);  % preallocate space for t and dimer, for speed xout = zeros(10,1e6); 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,1000) == 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';	

end

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

t = t(1:i); 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(7); 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);

function P = updateP(P0,xout,L)

for t = 1:L

   for i = 1:10

      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 f = graphP(P,v)

[m n] = size(P);

slab = {'P_tetR_Open' 'P_tetR_Closed' 'mRNA_tetR' 'TetR' 'TetR2' ...

   'P_cI_open' 'P_cI_closed' 'mRNA_cI' 'CI' 'CI2'};

i=0;

for j = [3 4 5 8 9 10]

   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)
   hold off

end