Code
From 2007.igem.org
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
