Team:Colombia/Scripting

From 2014.igem.org

(Difference between revisions)
(Created page with "{{Http://2014.igem.org/Team:Colombia Uniandes/Bootstrap}}")
 
(44 intermediate revisions not shown)
Line 1: Line 1:
-
{{Http://2014.igem.org/Team:Colombia Uniandes/Bootstrap}}
+
{{Http://2014.igem.org/Team:Colombia}}
 +
<html>
 +
<div class="span11" style="text-align: justify;">
 +
<br>
 +
<p>
 +
<br>
 +
<center><b> <font size="24">
 +
<br>
 +
Scripting
 +
</font><br>
 +
<h1  class="curs1"><font size="4">
 +
Feel free to expand and scroll through the text boxes in order to further examine the code.
 +
</font></h1></b></center>
 +
<br>
 +
<b><font color="#8A0808" size="3" >Deterministic Model</font></b>
 +
<br>
 +
This code creates the differential equations governing the concentration dinamics of each protein in our model, finds the steady state solutions and then solves them using the numerical aproximation method Runge-Kutta<br>
 +
<textarea readonly rows = "20">
 +
function y=det()
 +
%
 +
%--------------------------------------------%
 +
%                Runge Kutta 4 order aproximation          %
 +
%--------------------------------------------%
 +
%
 +
h=1200;                                                            % Maximum time
 +
m=0.005;                                                          % Step lenght [s]
 +
xi=[1,1,1,1,1,1,1,1,1];                                        %
 +
y=fsolve(@CondIni,xi);                                        %
 +
conInd=[1,1,1,1,1,1,1,1,0];                                %
 +
l=(0:m:h)';                                                          % Time vector
 +
%
 +
x=zeros(length(l),length(conInd));                        % Variable's Matrix, the variables in each column and time in each row
 +
C=zeros(1,length(l));                                            % zero vector to be filled
 +
x(1,:)=conInd;                                                      %
 +
%
 +
for k=1:length(l)-1
 +
    xk=x(k,:);                                                        % Capture of the last value in the matrix, the actual variable values
 +
    k1=difeq(l(k),xk);                                            % First slope of the RK4 method
 +
    k2=difeq(l(k)+m/2,xk+(m/2*k1)');                  % Second slope of the RK4 method
 +
    k3=difeq(l(k)+m/2,xk+(m/2*k2)');                  % Third slope of the RK4 method
 +
    k4=difeq(l(k)+m,xk+(m*k3)');                          % Fourth slope of the RK4 method
 +
    xk1=xk+m/6*(k1+2*k2+2*k3+k4)';                % New value calculation for the variables
 +
    %xk1=xk+m*difeq(l(k),xk)';                              % Another aproximation, Newton's method (comment line)                   
 +
    xk2=zeros(1,length(xk1));                                %
 +
%   
 +
    for p=1:length(xk1)
 +
        if(xk1(p)<0.00000001)
 +
            xk2(p)=0;
 +
        else   
 +
            xk2(p)=xk1(p);
 +
        end 
 +
    end
 +
    x(k+1,:)=xk2;                                                    % Variable vector actualization in the matrix
 +
end
 +
%
 +
for j=1:length(l)
 +
    if (l(j)>(300) & l(j)<(1000))
 +
        C(j)=2.25;
 +
    else
 +
        C(j)=0;
 +
    end
 +
end
 +
%
 +
 
 +
CS=x(:,1);                                                              %
 +
CSa=x(:,2);                                                            %
 +
Uf=x(:,3);                                                              %
 +
U=x(:,4);                                                              %
 +
Of=x(:,5);                                                              %
 +
O=x(:,6);                                                              %
 +
TR=x(:,7);                                                              %
 +
TA=x(:,8);                                                              %
 +
R=x(:,9);                                                                %
 +
%
 +
%
 +
figure(1)
 +
%
 +
hold on;
 +
plot(l,R,'r')
 +
plot(l,C)
 +
plot(l,Of,'g')
 +
%plot(l,TR)
 +
hold off;
 +
xlabel('Time')
 +
ylabel('Concetration (micromolar)')
 +
title('Activator response')
 +
%
 +
end
 +
%
 +
%
 +
function y=difeq(t,x)
 +
%
 +
global kcc kcd dsu duo fsu fuo gcs gcsa guf gu gof go gtr gta gr acs au ao ar atr ata btr bta br ho htr n
 +
%
 +
%--------------------------------------------%
 +
%                              VARIABLES                                %
 +
%--------------------------------------------%
 +
%
 +
C=input(t);                                    % Extracellular concentration of the cholerae autoinducer-1 (CAI-1)
 +
%
 +
CS=x(1);                                      % Concentration of the membrane bound CqsS protein (CAI-1 unbound=Inactive)
 +
CSa=x(2);                                    % Concentration of the membrane bound CqsS protein (CAI-1 bound=Active)
 +
Uf= x(3);                                      % Phosphorelay protein LuxU (phosphorylated) Concentration
 +
U=x(4);                                        % Phosphorelay protein LuxU (unphosphorylated) Concentration
 +
Of=x(5);                                      % Phosphorelay protein LuxO (phosphorylated) Concentration
 +
O=x(6);                                        % Phosphorelay protein LuxO (unphosphorylated) Concentration
 +
TR=x(7);                                      % Ptet Repressor protein concentration
 +
TA=x(8);                                      % Ptet Activator protein concentration
 +
R=x(9);                                        % Response molecule concentration
 +
%
 +
%--------------------------------------------%
 +
%                              PARAMETERS                            %
 +
%--------------------------------------------%
 +
%
 +
kcc=1;                                          % CAI1 and CqsS coupling Rate
 +
kcd=1;                                          % CAI1 and CqsS decoupling Rate
 +
dsu=4;                                          % LuxU* dephosphorylation rate through CqsS*
 +
duo=4;                                        % LuxO* dephosphorylation rate through LuxU
 +
fsu=2;                                          % LuxU* phosphorylation rate through CqsS
 +
fuo=2;                                          % LuxO* phosphorylation rate through LuxU*
 +
%
 +
gcs=1;                                          % CqsS protein decay rate
 +
gcsa=1;                                        % CqsS* protein decay rate
 +
guf=1;                                          % LuxU* protein decay rate
 +
gu=1;                                          % LuxU protein decay rate
 +
gof=1;                                          % LuxO* protein decay rate
 +
go=1;                                          % LuxO protein decay rate
 +
gtr=1;                                          % Ptet Repressor protein decay rate
 +
gta=1;                                          % Ptet Activator protein decay rate
 +
gr=1;                                          % Response molecule decay rate
 +
%
 +
acs=3;                                          % CS basal production rate
 +
au=3;                                          % LuxU basal production rate
 +
ao=3;                                          % LuxO basal production rate
 +
ar=0.01;                                      % response molecule basal production rate
 +
atr=0.01;                                      % TR basal production rate
 +
ata=0.01;                                      % TA basal production rate
 +
 
 +
btr=5;                                          % Maximum rate of TR expression
 +
bta=5;                                          % Maximum rate of TA expression
 +
br=5;                                          % Maximum rate of response molecule expression
 +
ho=1.5;                                      % LuxO*- DNA coupling rate
 +
htr=2;                                          % TRdomain-DNA coupling rate
 +
%
 +
n=1;                                            % Hill coefficient
 +
%
 +
%--------------------------------------------%
 +
%                                Equations                                %
 +
%--------------------------------------------%
 +
%
 +
dCS = acs + kcd*CSa - C*CS*kcc - gcs*CS;                                      % Differential equation governing the change in the concentration of the membrane bound CqsS protein (CAI-1 unbound=Inactive)through time
 +
dCSa = -kcd*CSa + C*CS*kcc - gcsa*CSa;                                        % Differential equation governing the change in the concentration of the membrane bound CqsS protein (CAI-1 bound=Active) through time
 +
dUf = U*Of*duo-Uf*O*fuo-CSa*Uf*dsu + CS*U*fsu - guf*Uf;          % Differential equation governing the change in the phosphorelay protein LuxU (phosphorylated) Concentration through time
 +
dU = au-U*Of*duo+Uf*O*fuo - CS*U*fsu + CSa*Uf*dsu -gu*U;      % Differential equation governing the change in the phosphorelay protein LuxU (unphosphorylated) Concentration through time
 +
dOf = Uf*O*fuo - U*Of*duo - gof*Of;                    % Differential equation governing the change in the phosphorelay protein LuxO (phosphorylated) Concentration through time
 +
dO = ao - Uf*O*fuo + U*Of*duo - go*O;                          % Differential equation governing the change in the phosphorelay protein LuxO (unphosphorylated) Concentration through time
 +
dTR = atr + (btr*Of^n)/(ho^n+Of^n) - gtr*TR;  % Differential equation governing the change in the ptet Repressor protein concentration through time
 +
dTA = ata + bta/((1+TR/TA)+(htr/TA)) - gta*TA; % Differential equation governing the change in the ptet Activator protein concentration through time
 +
dR = ar + br/((1+TR/TA)+(htr/TA)) - gr*R;      % Differential equation governing the change in the response molecule concentration through time
 +
%
 +
y(1)=dCS;                                      %
 +
y(2)=dCSa;                                    %
 +
y(3)=dUf;                                      %
 +
y(4)=dU;                                        %
 +
y(5)=dOf;                                      %
 +
y(6)=dO;                                        %
 +
y(7)=dTR;                                      %
 +
y(8)=dTA;                                      %
 +
y(9)=dR;                                        %
 +
%
 +
y=y';                                              % Transposing the vector because of matlab language restrictions
 +
%
 +
end
 +
%
 +
%
 +
function y=input(t)
 +
%--------------------------------------------%
 +
%                      CAI-1 (C) pulse simulation                %
 +
%--------------------------------------------%
 +
if t>50 && t<100 
 +
  y=10;
 +
elseif t>250 && t<1500
 +
    y=20;
 +
else
 +
    y=0;
 +
end
 +
%
 +
%
 +
</textarea>
 +
<br>
 +
<br>
 +
<b><font color="#8A0808" size="3" >Metropolis-Hastings Algorithm</font></b>
 +
<br>
 +
This code determines which set of values for missing parameters yield the most desirable response.<br>
 +
<textarea readonly rows = "20">
 +
function y=model(patr, pbtr, pho, pkttr, pata, pbta, phtr,pn)
 +
%
 +
global kcc kcd kcu kuc kuo kou kttr gcs gcsa guf gu gof go gtr gta gr gttr acs au ao ar atr ata bcu buc btr buo bou bta ho htr n
 +
%
 +
%
 +
%--------------------------------------------%
 +
%                              PARAMETERS                              %
 +
%--------------------------------------------%
 +
%
 +
kcc=1;                                                  % CAI1 and CqsS coupling Rate
 +
kcd=1;                                                  % CAI1 and CqsS decoupling Rate
 +
kcu=2;                                                  % Michaelis Menten Constant for the phoshporylation CqsS-LuxU
 +
kuc=1;                                                  % Michaelis Menten Constant of the phosphase process CqSS-LuxU
 +
kuo=2;                                                  % Michaelis Menten Constant for the phoshporylation LuxU-LuxO
 +
kou=1;                                                  % Michaelis Menten Constant of the phosphase process LuxU-LuxO
 +
kttr=pkttr;                                            %Dimerization rate of TetR
 +
%
 +
gcs=0.5;                                                % CqsS protein decay rate
 +
gcsa=1;                                                % CqsS* protein decay rate
 +
guf=0.12;                                              % LuxU* protein decay rate
 +
gu=0.65;                                              % LuxU protein decay rate
 +
gof=0.12;                                              % LuxO* protein decay rate
 +
go=0.12;                                              % LuxO protein decay rate
 +
gtr=0.12;                                              % Ptet Repressor protein decay rate
 +
gttr=0.12;                                            %Ptet2 dimer Repressor decay rate
 +
gta=0.12;                                              % Ptet Activator protein decay rate
 +
gr=0.12;                                                % Response molecule decay rate
 +
%
 +
acs=5;                                                  % CS basal production rate
 +
au=5;                                                    % LuxU basal production rate
 +
ao=5;                                                    % LuxO basal production rate
 +
ar=0;                                                    % response molecule basal production rate
 +
atr=patr;                                                % TR basal production rate
 +
ata=pata;                                              % TA basal production rate
 +
%
 +
bcu=1;                                                  % Phosphorylation rate of CqsS to LuxU
 +
buc=0.1;                                                % Phosphase rate of CqsS to LuxU
 +
buo=3.2;                                                % Phosphorylation rate LuxU-LuxO
 +
bou=1.58;                                              % Phosphase rate LuxU-LuxO
 +
btr=pbtr;                                                % Maximum rate of TR expression
 +
bta=pbta;                                              % Maximum rate of TA expression
 +
%
 +
ho=pho;                                                % LuxO*- DNA coupling rate
 +
htr=phtr;                                                % TRdomain-DNA coupling rate
 +
%
 +
n=pn;                                                    % Hill coefficient
 +
%
 +
%
 +
%
 +
%
 +
%
 +
%--------------------------------------------%
 +
%                Runge Kutta 4 order aproximation            %
 +
%--------------------------------------------%
 +
%
 +
h=210;                                                                                  % Maximum time
 +
m=1/12;                                                                              % Step lenght [h]
 +
conInd=[1,0,1,1,1,1,0,1,1,1];                                                %
 +
l=(0:m:h)';                                                                            % Time vector
 +
x=zeros(length(l),length(conInd));                                        % Variable's Matrix, the variables in each column and time in each row
 +
C=zeros(1,length(l));                                                            % zero vector to be filled
 +
x(1,:)=conInd;                                                                      %
 +
%
 +
for k=1:length(l)-1
 +
    xk=x(k,:);                                                                          % Capture of the last value in the matrix, the actual variable values
 +
    k1=difeq_Metropolis(l(k),xk);                                            % First slope of the RK4 method
 +
    k2=difeq_Metropolis(l(k)+m/2,xk+(m/2*k1)');                  % Second slope of the RK4 method
 +
    k3=difeq_Metropolis(l(k)+m/2,xk+(m/2*k2)');                  % Third slope of the RK4 method
 +
    k4=difeq_Metropolis(l(k)+m,xk+(m*k3)');                          % Fourth slope of the RK4 method
 +
    xk1=xk+m/6*(k1+2*k2+2*k3+k4)';                                  % New value calculation for the variables
 +
    %xk1=xk+m*difeq_Metropolis(l(k),xk)';                              % Another aproximation, Newton's method                     
 +
    xk2=zeros(1,length(xk1));                                                  %
 +
%   
 +
    for p=1:length(xk1)
 +
        if(xk1(p)<0.00000001)
 +
            xk2(p)=0;
 +
        else   
 +
            xk2(p)=xk1(p);
 +
        end 
 +
    end
 +
    x(k+1,:)=xk2;                                                  % Actualizacion del nuevo vector de variables en la matriz
 +
end
 +
%
 +
for j=1:length(l)
 +
    if (l(j)>100 && l(j)<200)
 +
        C(j)=20;
 +
    else
 +
        C(j)=0;
 +
           
 +
    end
 +
end
 +
%
 +
   
 +
CS=x(:,1);                                                            %
 +
CSa=x(:,2);                                                          %
 +
Uf=x(:,3);                                                            %
 +
U=x(:,4);                                                              %
 +
Of=x(:,5);                                                            %
 +
O=x(:,6);                                                              %
 +
TR=x(:,7);                                                            %                     
 +
TTR=x(:,8);                                                          %
 +
TA=x(:,9);                                                            %
 +
R=x(:,10);                                                            %
 +
 
 +
%figure;                                                                %
 +
%plot(l,C,l,R,'r');                                                    %
 +
%legend ('CAI','Chromoprotein');                            %
 +
%xlabel('Time [h]');                                                %
 +
%ylabel('Concentration [nanomolar]');                    %
 +
%title('Response');                                                  %
 +
%
 +
y=[R(95*12),R(100*12),R(105*12),R(110*12),R(115*12),R(120*12),R(125*12),R(130*12),R(135*12),R(140*12),R(145*12),R(150*12),R(155*12),R(160*12),R(165*12),R(170*12)];
 +
%
 +
%
 +
%
 +
%
 +
RFP=[3,3.1,3.15,3.5,4.8,6,11,15,16.5,18,18.5,19,19.5,20,20.25,20.5];
 +
t=[95,100,105,110,115,120,125,130,135,140,145,150,155,160,165,170];
 +
%
 +
steps=5000;
 +
%
 +
atr_walk=zeros(1,steps+1);
 +
btr_walk=zeros(1,steps+1);
 +
ho_walk=zeros(1,steps+1);
 +
kttr_walk=zeros(1,steps+1);
 +
ata_walk=zeros(1,steps+1);
 +
bta_walk=zeros(1,steps+1);
 +
htr_walk=zeros(1,steps+1);
 +
n_walk=zeros(1,steps+1);
 +
%
 +
atr_walk(1)=random('unif',0.5,1.5);
 +
btr_walk(1)=random('unif',4,6);
 +
ho_walk(1)=random('unif',1.2,1.8);
 +
kttr_walk(1)=random('unif',2,7);
 +
ata_walk(1)=random('unif',0.1,0.5);
 +
bta_walk(1)=random('unif',4,7);
 +
htr_walk(1)=random('unif',4,7);
 +
n_walk(1)=random('unif',1,3);
 +
%
 +
    for i=1:steps
 +
%
 +
    atr_prime=random('unif',atr_walk(i)-0.2,atr_walk(i)+0.2);
 +
    btr_prime=random('unif',btr_walk(i)-0.2,btr_walk(i)+0.2);
 +
    ho_prime=random('unif',ho_walk(i)-0.2,ho_walk(i)+0.2);
 +
    kttr_prime=random('unif',kttr_walk(i)-0.2,kttr_walk(i)+0.2);
 +
    ata_prime=random('unif',ata_walk(i)-0.2,ata_walk(i)+0.2);
 +
    bta_prime=random('unif',bta_walk(i)-0.2,bta_walk(i)+0.2);
 +
    htr_prime=random('unif',htr_walk(i)-0.2,htr_walk(i)+0.2);
 +
    n_prime=random('unif',n_walk(i)-0.2,n_walk(i)+0.2);
 +
% 
 +
    model_i=model(atr_walk(i),btr_walk(i),ho_walk(i),kttr_walk(i),ata_walk(i),bta_walk(i),htr_walk(i),n_walk(i));
 +
    model_prime=model(atr_prime,btr_prime,ho_prime,kttr_prime,ata_prime,bta_prime,htr_prime,n_prime);
 +
%
 +
    alpha=likelihood_lineal(RFP,model_prime)/likelihood_lineal(RFP,model_i);
 +
%   
 +
        if(alpha>=1)
 +
            atr_walk(i+1)=atr_prime;
 +
            btr_walk(i+1)=btr_prime;
 +
            ho_walk(i+1)=ho_prime;
 +
            kttr_walk(i+1)=kttr_prime;
 +
            ata_walk(i+1)=ata_prime;
 +
            bta_walk(i+1)=bta_prime;
 +
            htr_walk(i+1)=htr_prime;
 +
            n_walk(i+1)=n_prime;
 +
        else
 +
          % beta=random('unif',0,1);
 +
            %if(beta<=alpha)
 +
            %  atr_walk(i+1)=atr_prime;
 +
              % btr_walk(i+1)=btr_prime;
 +
              %ho_walk(i+1)=ho_prime;
 +
              %kttr_walk(i+1)=kttr_prime;
 +
              %ata_walk(i+1)=ata_prime;
 +
              %bta_walk(i+1)=bta_prime;
 +
              %htr_walk(i+1)=htr_prime;
 +
              %n_walk(i+1)=n_prime;
 +
            %else
 +
              atr_walk(i+1)=atr_walk(i);
 +
              btr_walk(i+1)=btr_walk(i);
 +
              ho_walk(i+1)=ho_walk(i);
 +
              kttr_walk(i+1)=kttr_walk(i);
 +
              ata_walk(i+1)=ata_walk(i);
 +
              bta_walk(i+1)=bta_walk(i);
 +
              htr_walk(i+1)=htr_walk(i);
 +
              n_walk(i+1)=n_walk(i);
 +
            %end
 +
        end
 +
    i
 +
    end
 +
% 
 +
figure (1)
 +
plot(t,RFP,'b',t,model(atr_walk(steps+1),btr_walk(steps+1),ho_walk(steps+1),kttr_walk(steps+1),ata_walk(steps+1),bta_walk(steps+1),htr_walk(steps+1),n_walk(steps+1)),'r')
 +
xlabel('tiempo')
 +
ylabel('Reportero')
 +
title('Respuesta')
 +
%
 +
figure (2)
 +
scatter((1:steps+1),atr_walk)
 +
xlabel('step')
 +
ylabel('valor')
 +
title('atr')
 +
%
 +
figure (3)
 +
scatter((1:steps+1),btr_walk)
 +
xlabel('step')
 +
ylabel('valor')
 +
title('btr')
 +
%
 +
figure (4)
 +
scatter((1:steps+1),ho_walk)
 +
xlabel('step')
 +
ylabel('valor')
 +
title('ho')
 +
%
 +
figure (5)
 +
scatter((1:steps+1),kttr_walk)
 +
xlabel('step')
 +
ylabel('valor')
 +
title('kttr')
 +
%
 +
figure (6)
 +
scatter((1:steps+1),ata_walk)
 +
xlabel('step')
 +
ylabel('valor')
 +
title('ata')
 +
%
 +
figure (7)
 +
scatter((1:steps+1),bta_walk)
 +
xlabel('step')
 +
ylabel('valor')
 +
title('bta')
 +
%
 +
figure (8)
 +
scatter((1:steps+1),htr_walk)
 +
xlabel('step')
 +
ylabel('valor')
 +
title('htr')
 +
%
 +
figure (9)
 +
scatter((1:steps+1),n_walk)
 +
xlabel('step')
 +
ylabel('valor')
 +
title('n')
 +
%
 +
%
 +
function y=input_Metropolis(t)
 +
%--------------------------------------------%
 +
%                    CAI-1 (C) pulse simulation                    %
 +
%--------------------------------------------%
 +
if t>100 && t<200 
 +
  y=20;
 +
else
 +
    y=0;
 +
end
 +
</textarea>
 +
<br>
 +
<br>
 +
<b><font color="#8A0808" size="3" >Sensitivity Analysis</font></b>
 +
<br>
 +
This code points out which are the critical parameters in the system (those that change the response drastically).<br>
 +
<textarea readonly rows = "20">
 +
%
 +
global kcc kcd kcu kuc kuo kou kttr gcs gcsa guf gu gof go gtr gta gr gttr acs au ao ar atr ata bcu buc btr buo bou bta ho htr n
 +
%
 +
%
 +
%--------------------------------------------%
 +
%                                PARAMETERS                            %
 +
%--------------------------------------------%
 +
%
 +
kcc=1;                                                          % CAI1 and CqsS coupling Rate
 +
kcd=1;                                                          % CAI1 and CqsS decoupling Rate
 +
kcu=2;                                                          % Michaelis Menten Constant for the phoshporylation CqsS-LuxU
 +
kuc=1;                                                          % Michaelis Menten Constant of the phosphase process CqSS-LuxU
 +
kuo=2;                                                          % Michaelis Menten Constant for the phoshporylation LuxU-LuxO
 +
kou=1;                                                          % Michaelis Menten Constant of the phosphase process LuxU-LuxO
 +
kttr=5;                                                          % Dimerization rate of TetR ---[Sensitivity Analysis]--- [5-50/5]
 +
%
 +
gcs=0.5;                                                        % CqsS protein decay rate
 +
gcsa=1;                                                          % CqsS* protein decay rate
 +
guf=0.12;                                                      % LuxU* protein decay rate
 +
gu=0.65;                                                        % LuxU protein decay rate
 +
gof=0.12;                                                      % LuxO* protein decay rate
 +
go=0.12;                                                        % LuxO protein decay rate
 +
gtr=0.12;                                                      % Ptet Repressor protein decay rate
 +
gttr=0.12;                                                      % Ptet2 dimer Repressor decay rate
 +
gta=0.12;                                                      % Ptet Activator protein decay rate
 +
gr=0.12;                                                        % Response molecule decay rate
 +
%
 +
acs=5;                                                            % CS basal production rate
 +
au=5;                                                              % LuxU basal production rate
 +
ao=5;                                                              % LuxO basal production rate
 +
ar=0;                                                              % response molecule basal production rate
 +
atr=1;                                                              % TR basal production rate  ---[Sensitivity Analysis]--- [0.5-5/0.5]
 +
ata=0.1;                                                          % TA basal production rate ---[Sensitivity Analysis]--- [0.1-10/0.5]
 +
%
 +
bcu=1;                                                            % Phosphorylation rate of CqsS to LuxU
 +
buc=0.1;                                                          % Phosphase rate of CqsS to LuxU
 +
buo=3.2;                                                          % Phosphorylation rate LuxU-LuxO
 +
bou=1.58;                                                        % Phosphase rate LuxU-LuxO
 +
btr=5;                                                              % Maximum rate of TR expression ---[Sensitivity Analysis]--- [1-10/1]
 +
bta=5;                                                              % Maximum rate of TA expression ---[Sensitivity Analysis]--- [1-10/1]
 +
%
 +
ho=1.5;                                                            % LuxO*- DNA coupling rate ---[Sensitivity Analysis]--- [0.5-2.5/0.5]
 +
htr=5;                                                              % TRdomain-DNA coupling rate ---[Sensitivity Analysis]--- [1-10/1]
 +
%
 +
n=2;                                                                  % Hill's coefficient ---[Sensitivity Analysis]--- [1-4/0.5]
 +
%
 +
h=400;                                                              % Maximum time
 +
m=1/10;                                                            % Step lenght [h]
 +
%
 +
param= (5:50:5);                                                % Range for the parameter to vary, evaluate at least 2 orders of magnitude chose maximum 20 points.
 +
rta=zeros(3,length(param));                                % Array that saves the canche in concentration of LuxOF, the response and time of response.
 +
for w=1:length(param)                                        % Varying the value of the parameter
 +
    kttr= param(w);                                              % New parameter value (This MUST be changed everytime a parameter is being evaluated)
 +
%
 +
%
 +
%--------------------------------------------%
 +
%              Runge Kutta 4 order aproximation            %
 +
%--------------------------------------------%
 +
%
 +
%
 +
%
 +
conInd=[1,0,1,1,1,1,0,1,1,1];                              %
 +
l=(0:m:h)';                                                            % Time vector
 +
%
 +
x=zeros(length(l),length(conInd));                        % Variable's Matrix, the variables in each column and time in each row
 +
C=zeros(1,length(l));                                            % zero vector to be filled
 +
x(1,:)=conInd;                                                      %
 +
%
 +
for k=1:length(l)-1
 +
    xk=x(k,:);                                                          % Capture of the last value in the matrix, the actual variable values
 +
    k1=difeq(l(k),xk);                                              % First slope of the RK4 method
 +
    k2=difeq(l(k)+m/2,xk+(m/2*k1)');                    % Second slope of the RK4 method
 +
    k3=difeq(l(k)+m/2,xk+(m/2*k2)');                    % Third slope of the RK4 method
 +
    k4=difeq(l(k)+m,xk+(m*k3)');                            % Fourth slope of the RK4 method
 +
    xk1=xk+m/6*(k1+2*k2+2*k3+k4)';                  % New value calculation for the variables
 +
    xk2=zeros(1,length(xk1));                                  %
 +
%   
 +
    for p=1:length(xk1)
 +
        if(xk1(p)<0.00000001)
 +
            xk2(p)=0;
 +
        else   
 +
            xk2(p)=xk1(p);
 +
        end 
 +
    end
 +
    x(k+1,:)=xk2;                                % Actualizacion del nuevo vector de variables en la matriz
 +
end
 +
%
 +
R=x(:,10);
 +
Of=x(:,5);%
 +
%
 +
%
 +
it=141;
 +
%
 +
  for it=1011:3001
 +
%
 +
%
 +
  if (R(it)-R(1011))/(R(1011))> 1.5      % If the response changes 1.5-fold relative to previous iteration
 +
%
 +
  rta(1,w)= l(it)-100;
 +
%
 +
%
 +
%
 +
  end
 +
%
 +
%
 +
%
 +
  end 
 +
%
 +
%
 +
%
 +
%
 +
rta(2,w)= Of(2701)-Of(801);
 +
rta(3,w)= R(2701)-R(801);
 +
%
 +
end
 +
%
 +
%Here Labels should be changed for them to make proper sense for each evaluated parameter
 +
%
 +
figure (1)                             
 +
plot(param,rta(1,:),'r')
 +
xlabel('Dimerization rate of TetR');                                   
 +
ylabel('Time of response')
 +
title('Sensitivity analysis')
 +
%
 +
figure (2)
 +
plot(param,rta(2,:),'g')
 +
xlabel('Dimerization rate of TetR');
 +
ylabel('Change in LuxO Phosphorylated')
 +
title('Sensitivity analysis')
 +
%
 +
figure (3)
 +
plot(param,rta(3,:)) 
 +
xlabel('Dimerization rate of TetR');
 +
ylabel('Change in response (nM)')
 +
title('Sensitivity analysis')
 +
</textarea>
 +
<br>
 +
<br>
 +
<b><font color="#8A0808" size="3" >Stochastic Model</font></b>
 +
<br>
 +
Sometimes probabilistic models better describe certain systems, This required a fairly big amount of computational power<br>
 +
<textarea readonly rows = "20">
 +
%Stochastic Simulation for the system
 +
%
 +
%%----Parameters of the System----
 +
global kcc kcd kcu kuc kuo kou kttr gcs gcsa guf gu gof go gtr gta gr acs au ao ar atr ata bcu buc btr buo bou bta ho htr n gttr
 +
%
 +
%
 +
numCells = 20;
 +
numEvents= 50000;
 +
%
 +
%
 +
%
 +
%
 +
kcc=1;                                              % CAI1 and CqsS coupling Rate
 +
kcd=1;                                                % CAI1 and CqsS decoupling Rate
 +
kcu=1.2044000;                                        % Michaelis Menten Constant for the phoshporylation CqsS-LuxU
 +
kuc=0.6022000;                                          % Michaelis Menten Constant of the phosphase process CqSS-LuxU
 +
kuo=1.2044000;                                        % Michaelis Menten Constant for the phoshporylation LuxU-LuxO
 +
kou=0.6022000;                                        % Michaelis Menten Constant of the phosphase process LuxU-LuxO
 +
kttr=5;                                              %Dimerization rate of TetR
 +
%
 +
gcs=0.500;                                              % CqsS protein decay rate
 +
gcsa=1;                                            % CqsS* protein decay rate
 +
guf=0.120;                                              % LuxU* protein decay rate
 +
gu=0.650;                                              % LuxU protein decay rate
 +
gof=0.120;                                              % LuxO* protein decay rate
 +
go=0.650;                                              % LuxO protein decay rate
 +
gtr=0.120;                                              % Ptet Repressor protein decay rate
 +
gttr=0.120;                                              %Ptet2 dimer Repressor decay rate
 +
gta=0.120;                                              % Ptet Activator protein decay rate
 +
gr=0.050;                                              % Response molecule decay rate
 +
%
 +
acs=5;                                              % CS basal production rate
 +
au=5;                                              % LuxU basal production rate
 +
ao=5;                                              % LuxO basal production rate
 +
ar=0.1;                                            % response molecule basal production rate
 +
atr=1;                                          % TR basal production rate
 +
ata=0.1;                                          % TA basal production rate
 +
%
 +
bcu=1;                                              % Phosphorylation rate of CqsS to LuxU
 +
buc=0.1;                                              % Phosphase rate of CqsS to LuxU
 +
buo=3.2;                                              % Phosphorylation rate LuxU-LuxO
 +
bou=1.58;                                              % Phosphase rate LuxU-LuxO
 +
btr=5;                                              % Maximum rate of TR expression
 +
bta=5;                                              % Maximum rate of TA expression
 +
%
 +
ho=1.5;                                              % LuxO*- DNA coupling rate
 +
htr=5;                                              % TRdomain-DNA coupling rate
 +
%
 +
n=3;% Hill coefficient
 +
%
 +
%%----Stochastic Simulation----
 +
%
 +
%Variables of the System
 +
C=zeros(numCells,numEvents);
 +
CS=zeros(numCells,numEvents);
 +
CSa=zeros(numCells,numEvents);
 +
Uf=zeros(numCells,numEvents);
 +
U=zeros(numCells,numEvents);
 +
Of=zeros(numCells,numEvents);
 +
O=zeros(numCells,numEvents);
 +
TR=zeros(numCells,numEvents);
 +
TRR=zeros(numCells,numEvents);
 +
TA=zeros(numCells,numEvents);
 +
R=zeros(numCells,numEvents);
 +
%
 +
C(:,1)=0;
 +
CS(:,1)=10;
 +
CSa(:,1)=1;
 +
Uf(:,1)=15;
 +
U(:,1)=5;
 +
Of(:,1)=40;
 +
O(:,1)=5;
 +
TA(:,1)=10;
 +
TR(:,1)=2;
 +
TRR(:,1)=1;
 +
%
 +
%
 +
%
 +
%
 +
t = zeros(numCells,numEvents); %Time
 +
tr = zeros(numCells,numEvents);
 +
%
 +
%
 +
%Stochastic Simulation
 +
%
 +
for i =1:numCells
 +
   
 +
    i
 +
    for j = 1:numEvents
 +
       
 +
        j
 +
       
 +
        %Random numbers drawn from uniform distribution between zero and
 +
        %one
 +
       
 +
        rand1 = rand(1); %This determine time
 +
        rand2 = rand(1); %This determine Event
 +
       
 +
        %Events (Description in attached document)
 +
       
 +
        e1= kcc*CS(i,j)*C(i,j);
 +
        e2= kcd*CSa(i,j);
 +
        e3= gcs*CS(i,j);
 +
        e4= acs;
 +
        e5= gcsa*CSa(i,j);
 +
        e6= au;
 +
        e7= (bcu*U(i,j)*CS(i,j))/(kcu+U(i,j));
 +
        e8= (buc*Uf(i,j)*CSa(i,j))/(kuc+Uf(i,j));
 +
        e9= (buo*O(i,j)*Uf(i,j))/(kuo+O(i,j));
 +
        e10=(bou*Of(i,j)*U(i,j))/(kou+Of(i,j));
 +
        e11=gu*U(i,j);
 +
        e12=guf*Uf(i,j);
 +
        e13=ao;
 +
        e14=go*O(i,j);
 +
        e15=gof*Of(i,j);
 +
       
 +
        e16=atr;
 +
        e17=(btr*Of(i,j)^n)/((ho.^n)+(Of(i,j).^n));
 +
        e18=kttr*TR(i,j)^2;
 +
        e19=gtr*TR(i,j);
 +
       
 +
        e20=ata;
 +
        e21=ar;
 +
        e22=bta/(1+(TRR(i,j)/TA(i,j))+(htr/TA(i,j)));
 +
        e23=gta*TA(i,j);
 +
        e24=gr*R(i,j);
 +
        e25=TRR(i,j)*gttr;
 +
       
 +
       
 +
        %Vector of Events
 +
        vecEvents = [e1,e2,e3,e4,e5,e6,e7,e8,e9,e10,e11,e12,e13,e14,e15,e16,e17,e18,e19,e20,e21,e22,e23,e24,e25];
 +
       
 +
        %Sum of Events
 +
        sumEvents = sum(vecEvents);
 +
       
 +
        %Waiting time until the next ocurrence of a reaction time
 +
        %Probabilistic Distribution
 +
        tau= (1/sumEvents)*log(1/rand1);
 +
        tr(i,j)=tau;
 +
        t(i,j+1)= t(i,j)+ tau;
 +
       
 +
        %CAI outside the cell
 +
       
 +
        C(i,j+1)=funcImpulso(t(i,j+1));
 +
       
 +
        %if(t(i,j)<130)
 +
        %  C(i,j)=funcImpulso(t(i,j+1))*957/2575;
 +
        %end
 +
       
 +
        %Simulation Creation/Destruction of each molecule
 +
       
 +
        if (rand2<(vecEvents(1)/sumEvents))  %CqsS activation
 +
           
 +
            if (CS(i,j)<=0)
 +
               
 +
                CS(i,j+1)=0;
 +
            else
 +
                CS(i,j+1)= CS(i,j)-1;
 +
            end
 +
           
 +
            CSa(i,j+1)= CSa(i,j)+1;
 +
            U(i,j+1)= U(i,j);
 +
            Uf(i,j+1)= Uf(i,j);
 +
            O(i,j+1)= O(i,j);
 +
            Of(i,j+1)= Of(i,j);
 +
            TR(i,j+1)= TR(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TA(i,j+1)= TA(i,j);
 +
            R(i,j+1)= R(i,j);
 +
           
 +
           
 +
           
 +
        elseif (rand2< (sum(vecEvents(1:2)))/sumEvents) %CqsS desactivation
 +
           
 +
            if (CSa(i,j)<=0)
 +
               
 +
                CSa(i,j+1)= 0;
 +
               
 +
            else
 +
                CSa(i,j+1)= CSa(i,j)-1;
 +
            end
 +
           
 +
            CS(i,j+1)= CS(i,j)+1;
 +
            U(i,j+1)= U(i,j);
 +
            Uf(i,j+1)= Uf(i,j);
 +
            O(i,j+1)= O(i,j);
 +
            Of(i,j+1)= Of(i,j);
 +
            TR(i,j+1)= TR(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TA(i,j+1)= TA(i,j);
 +
            R(i,j+1)= R(i,j);
 +
           
 +
           
 +
           
 +
        elseif (rand2< (sum(vecEvents(1:3)))/sumEvents) %CqsS degradation
 +
           
 +
            if (CS(i,j)<=0)
 +
               
 +
                CS(i,j+1)=0;
 +
            else
 +
               
 +
                CS(i,j+1)= CS(i,j)-1;
 +
            end
 +
           
 +
            CSa(i,j+1)= CSa(i,j);
 +
            U(i,j+1)= U(i,j);
 +
            Uf(i,j+1)= Uf(i,j);
 +
            O(i,j+1)= O(i,j);
 +
            Of(i,j+1)= Of(i,j);
 +
            TR(i,j+1)= TR(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TA(i,j+1)= TA(i,j);
 +
            R(i,j+1)= R(i,j);
 +
           
 +
        elseif (rand2< (sum(vecEvents(1:4)))/sumEvents)  %CqsS basal production
 +
           
 +
           
 +
            CS(i,j+1)= CS(i,j)+1;
 +
            CSa(i,j+1)= CSa(i,j);
 +
            U(i,j+1)= U(i,j);
 +
            Uf(i,j+1)= Uf(i,j);
 +
            O(i,j+1)= O(i,j);
 +
            Of(i,j+1)= Of(i,j);
 +
            TR(i,j+1)= TR(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TA(i,j+1)= TA(i,j);
 +
            R(i,j+1)= R(i,j);
 +
           
 +
        elseif (rand2< (sum(vecEvents(1:5)))/sumEvents) %CqsS* degradation
 +
           
 +
            if (CSa(i,j)<=0)
 +
               
 +
                CSa(i,j+1)=0;
 +
            else
 +
               
 +
                CSa(i,j+1)= CSa(i,j)-1;
 +
            end
 +
           
 +
            CS(i,j+1)= CS(i,j);
 +
            U(i,j+1)= U(i,j);
 +
            Uf(i,j+1)= Uf(i,j);
 +
            O(i,j+1)= O(i,j);
 +
            Of(i,j+1)= Of(i,j);
 +
            TR(i,j+1)= TR(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TA(i,j+1)= TA(i,j);
 +
            R(i,j+1)= R(i,j);
 +
           
 +
           
 +
        elseif (rand2< (sum(vecEvents(1:6)))/sumEvents)  %LuxU basal production
 +
           
 +
           
 +
            U(i,j+1)= U(i,j)+1;
 +
            CS(i,j+1)= CS(i,j);
 +
            CSa(i,j+1)= CSa(i,j);
 +
            Uf(i,j+1)= Uf(i,j);
 +
            O(i,j+1)= O(i,j);
 +
            Of(i,j+1)= Of(i,j);
 +
            TR(i,j+1)= TR(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TA(i,j+1)= TA(i,j);
 +
            R(i,j+1)= R(i,j);
 +
           
 +
           
 +
        elseif (rand2< (sum(vecEvents(1:7)))/sumEvents)  %LuxU fosforilation
 +
           
 +
            if (U(i,j)<=0)
 +
               
 +
                U(i,j+1)=0;
 +
            else
 +
                U(i, j+1)= U(i,j)-1;
 +
            end
 +
           
 +
            CS(i,j+1)= CS(i,j);
 +
            CSa(i,j+1)= CSa(i,j);
 +
            Uf(i,j+1)= Uf(i,j)+1;
 +
            O(i,j+1)= O(i,j);
 +
            Of(i,j+1)= Of(i,j);
 +
            TR(i,j+1)= TR(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TA(i,j+1)= TA(i,j);
 +
            R(i,j+1)= R(i,j);
 +
           
 +
           
 +
        elseif (rand2< (sum(vecEvents(1:8)))/sumEvents)  %LuxU dephosphorilation
 +
           
 +
            if (Uf(i,j)<=0)
 +
               
 +
                Uf(i,j+1)=0;
 +
            else
 +
                Uf(i, j+1)= Uf(i,j)-1;
 +
            end
 +
           
 +
            CS(i,j+1)= CS(i,j);
 +
            CSa(i,j+1)= CSa(i,j);
 +
            U(i,j+1)= U(i,j)+1;
 +
            O(i,j+1)= O(i,j);
 +
            Of(i,j+1)= Of(i,j);
 +
            TR(i,j+1)= TR(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TA(i,j+1)= TA(i,j);
 +
            R(i,j+1)= R(i,j);
 +
           
 +
           
 +
        elseif (rand2< (sum(vecEvents(1:9)))/sumEvents)  %LuxO phosphorilation
 +
           
 +
            if (Uf(i,j)<=0 || O(i,j)<=0 )
 +
                if(Uf(i,j)<0)
 +
                    Uf(i,j+1)=0;
 +
                else
 +
                    Uf(i,j+1)=Uf(i,j);
 +
                end
 +
                if(O(i,j)<0)
 +
                    O(i,j+1)=0;
 +
                else
 +
                    O(i,j+1)=O(i,j);
 +
                end
 +
             
 +
            else
 +
                Uf(i, j+1)= Uf(i,j)-1;
 +
                O(i,j+1)=O(i,j)-1;
 +
            end
 +
           
 +
            CS(i,j+1)= CS(i,j);
 +
            CSa(i,j+1)= CSa(i,j);
 +
            U(i,j+1)= U(i,j)+1;
 +
            Of(i,j+1)= Of(i,j)+1;
 +
            TR(i,j+1)= TR(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TA(i,j+1)= TA(i,j);
 +
            R(i,j+1)= R(i,j);
 +
           
 +
        elseif (rand2< (sum(vecEvents(1:10)))/sumEvents)  %LuxO dephosphorilation
 +
           
 +
            if (U(i,j)<=0)
 +
               
 +
                U(i,j+1)=0;
 +
            else
 +
                U(i,j+1)= U(i,j)-1;
 +
            end
 +
           
 +
            if (Of(i,j)<=0)
 +
               
 +
                Of(i,j+1)=0;
 +
            else
 +
                Of(i, j+1)= Of(i,j)-1;
 +
            end
 +
           
 +
           
 +
            CS(i,j+1)= CS(i,j);
 +
            CSa(i,j+1)= CSa(i,j);
 +
            Uf(i,j+1)= Uf(i,j)+1;
 +
            O(i,j+1)= O(i,j)+1;
 +
            TR(i,j+1)= TR(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TA(i,j+1)= TA(i,j);
 +
            R(i,j+1)= R(i,j);
 +
           
 +
        elseif (rand2< (sum(vecEvents(1:11)))/sumEvents)  %LuxU degradation
 +
           
 +
            if (U(i,j)<=0)
 +
               
 +
                U(i,j+1)=0;
 +
            else
 +
                U(i, j+1)= U(i,j)-1;
 +
            end
 +
           
 +
            CS(i,j+1)= CS(i,j);
 +
            CSa(i,j+1)= CSa(i,j);
 +
            Uf(i,j+1)= Uf(i,j);
 +
            O(i,j+1)= O(i,j);
 +
            Of(i,j+1)= Of(i,j);
 +
            TR(i,j+1)= TR(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TA(i,j+1)= TA(i,j);
 +
            R(i,j+1)= R(i,j);
 +
     
 +
        elseif (rand2< (sum(vecEvents(1:12)))/sumEvents )  %LuxU* degradation
 +
           
 +
            if (Uf(i,j)<=0)
 +
               
 +
                Uf(i,j+1)=0;
 +
            else
 +
                Uf(i, j+1)= Uf(i,j)-1;
 +
            end
 +
           
 +
            CS(i,j+1)= CS(i,j);
 +
            CSa(i,j+1)= CSa(i,j);
 +
            U(i,j+1)= U(i,j);
 +
            O(i,j+1)= O(i,j);
 +
            Of(i,j+1)= Of(i,j);
 +
            TR(i,j+1)= TR(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TA(i,j+1)= TA(i,j);
 +
            R(i,j+1)= R(i,j);
 +
           
 +
           
 +
        elseif (rand2< (sum(vecEvents(1:13)))/sumEvents )%Basal production LuxO
 +
           
 +
           
 +
            U(i, j+1)= U(i,j);
 +
            CS(i,j+1)= CS(i,j);
 +
            CSa(i,j+1)= CSa(i,j);
 +
            Uf(i,j+1)= Uf(i,j);
 +
            O(i,j+1)= O(i,j)+1;
 +
            Of(i,j+1)= Of(i,j);
 +
            TR(i,j+1)= TR(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TA(i,j+1)= TA(i,j);
 +
            R(i,j+1)= R(i,j);
 +
           
 +
           
 +
        elseif (rand2< (sum(vecEvents(1:14)))/sumEvents) %LuxO degradation
 +
           
 +
            if (O(i,j)<=0)
 +
               
 +
                O(i,j+1)=0;
 +
            else
 +
                O(i, j+1)= O(i,j)-1;
 +
            end
 +
           
 +
            CS(i,j+1)= CS(i,j);
 +
            CSa(i,j+1)= CSa(i,j);
 +
            U(i,j+1)= U(i,j);
 +
            Uf(i,j+1)= Uf(i,j);
 +
            Of(i,j+1)= Of(i,j);
 +
            TR(i,j+1)= TR(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TA(i,j+1)= TA(i,j);
 +
            R(i,j+1)= R(i,j);
 +
           
 +
        elseif (rand2< (sum(vecEvents(1:15)))/sumEvents ) %LuxO* degradation
 +
           
 +
            if (Of(i,j)<=0)
 +
               
 +
                Of(i,j+1)=0;
 +
            else
 +
                Of(i, j+1)= Of(i,j)-1;
 +
            end
 +
           
 +
            CS(i,j+1)= CS(i,j);
 +
            CSa(i,j+1)= CSa(i,j);
 +
            U(i,j+1)= U(i,j);
 +
            O(i,j+1)= O(i,j);
 +
            Uf(i,j+1)= Uf(i,j);
 +
            TR(i,j+1)= TR(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TA(i,j+1)= TA(i,j);
 +
            R(i,j+1)= R(i,j);
 +
           
 +
        elseif (rand2< (sum(vecEvents(1:16)))/sumEvents) %Basal Production tetR
 +
           
 +
           
 +
            TR(i, j+1)= TR(i,j)+1;
 +
            CS(i,j+1)= CS(i,j);
 +
            CSa(i,j+1)= CSa(i,j);
 +
            U(i,j+1)= U(i,j);
 +
            O(i,j+1)= O(i,j);
 +
            Of(i,j+1)= Of(i,j);
 +
            Uf(i,j+1)= Uf(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TA(i,j+1)= TA(i,j);
 +
            R(i,j+1)= R(i,j);
 +
           
 +
        elseif (rand2< (sum(vecEvents(1:17)))/sumEvents )%Active production tetR
 +
           
 +
           
 +
            TR(i, j+1)= TR(i,j)+1;
 +
            CS(i,j+1)= CS(i,j);
 +
            CSa(i,j+1)= CSa(i,j);
 +
            U(i,j+1)= U(i,j);
 +
            O(i,j+1)= O(i,j);
 +
            Of(i,j+1)= Of(i,j);
 +
            Uf(i,j+1)= Uf(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TA(i,j+1)= TA(i,j);
 +
            R(i,j+1)= R(i,j);
 +
         
 +
           
 +
        elseif (rand2< (sum(vecEvents(1:18)))/sumEvents) %tetR dimerization 
 +
           
 +
            if (TR(i,j)<=0)
 +
               
 +
                TR(i,j+1)=0;
 +
                TRR(i, j+1)= TRR(i,j);
 +
               
 +
            elseif (TR(i,j)==1)
 +
                TR(i, j+1)=1;
 +
                TRR(i, j+1)= TRR(i,j);
 +
            else
 +
                TR(i, j+1)= TR(i,j)-2;
 +
                TRR(i, j+1)= TRR(i,j)+1;
 +
            end
 +
           
 +
            CS(i,j+1)= CS(i,j);
 +
            CSa(i,j+1)= CSa(i,j);
 +
            U(i,j+1)= U(i,j);
 +
            O(i,j+1)= O(i,j);
 +
            Of(i,j+1)= Of(i,j);
 +
            Uf(i,j+1)= Uf(i,j);
 +
            TA(i,j+1)= TA(i,j);
 +
            R(i,j+1)= R(i,j);
 +
           
 +
        elseif (rand2< (sum(vecEvents(1:19)))/sumEvents) %tetR degradation 
 +
           
 +
            if (TR(i,j)<=0)
 +
               
 +
                TR(i,j+1)=0;
 +
            else
 +
                TR(i,j+1)= TR(i,j)-1;
 +
            end
 +
           
 +
            CS(i,j+1)= CS(i,j);
 +
            CSa(i,j+1)= CSa(i,j);
 +
            U(i,j+1)= U(i,j);
 +
            O(i,j+1)= O(i,j);
 +
            Of(i,j+1)= Of(i,j);
 +
            Uf(i,j+1)= Uf(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TA(i,j+1)= TA(i,j);
 +
            R(i,j+1)= R(i,j);
 +
       
 +
        elseif (rand2< (sum(vecEvents(1:20)))/sumEvents )%Basal production Tet activator 
 +
           
 +
           
 +
            TA(i, j+1)= TA(i,j)+1;
 +
            CS(i,j+1)= CS(i,j);
 +
            CSa(i,j+1)= CSa(i,j);
 +
            U(i,j+1)= U(i,j);
 +
            O(i,j+1)= O(i,j);
 +
            Of(i,j+1)= Of(i,j);
 +
            Uf(i,j+1)= Uf(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TR(i,j+1)= TR(i,j);
 +
            R(i,j+1)= R(i,j);
 +
       
 +
        elseif (rand2< (sum(vecEvents(1:21)))/sumEvents) %Basal production Reporter protein 
 +
           
 +
 
 +
            R(i, j+1)= R(i,j)+1;
 +
            CS(i,j+1)= CS(i,j);
 +
            CSa(i,j+1)= CSa(i,j);
 +
            U(i,j+1)= U(i,j);
 +
            O(i,j+1)= O(i,j);
 +
            Of(i,j+1)= Of(i,j);
 +
            Uf(i,j+1)= Uf(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TR(i,j+1)= TR(i,j);
 +
            TA(i,j+1)= TA(i,j);
 +
           
 +
        elseif (rand2< (sum(vecEvents(1:22)))/sumEvents )%Active production Reporter protein & activator 
 +
           
 +
            TA(i,j+1)= TA(i,j)+1;
 +
            R(i,j+1)= R(i,j)+1;
 +
            CS(i,j+1)= CS(i,j);
 +
            CSa(i,j+1)= CSa(i,j);
 +
            U(i,j+1)= U(i,j);
 +
            O(i,j+1)= O(i,j);
 +
            Of(i,j+1)= Of(i,j);
 +
            Uf(i,j+1)= Uf(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TR(i,j+1)= TR(i,j);
 +
           
 +
        elseif (rand2< (sum(vecEvents(1:23)))/sumEvents) %Tet activator degradation 
 +
           
 +
            if (TA(i,j)<=0)
 +
               
 +
                TA(i,j+1)=0;
 +
            else
 +
                TA(i, j+1)= TA(i,j)-1;
 +
            end
 +
           
 +
            CS(i,j+1)= CS(i,j);
 +
            CSa(i,j+1)= CSa(i,j);
 +
            U(i,j+1)= U(i,j);
 +
            O(i,j+1)= O(i,j);
 +
            Of(i,j+1)= Of(i,j);
 +
            Uf(i,j+1)= Uf(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TR(i,j+1)= TR(i,j);
 +
            R(i,j+1)= R(i,j);
 +
           
 +
        elseif (rand2 < (sum(vecEvents(1:24)))/sumEvents) %Reporter protein activator degradation 
 +
           
 +
            if (R(i,j)<=0)
 +
                R(i,j+1)=0;
 +
            else
 +
                R(i, j+1)= R(i,j)-1;
 +
            end
 +
           
 +
            CS(i,j+1)= CS(i,j);
 +
            CSa(i,j+1)= CSa(i,j);
 +
            U(i,j+1)= U(i,j);
 +
            O(i,j+1)= O(i,j);
 +
            Of(i,j+1)= Of(i,j);
 +
            Uf(i,j+1)= Uf(i,j);
 +
            TRR(i,j+1)= TRR(i,j);
 +
            TR(i,j+1)= TR(i,j);
 +
            TA(i,j+1)= TA(i,j);
 +
           
 +
       
 +
        elseif (rand2 < (sum(vecEvents(1:25)))/sumEvents) %TetR dimers degradation 
 +
           
 +
            if (TRR(i,j)<=0)
 +
                TRR(i,j+1)=0;
 +
            else
 +
                TRR(i, j+1)= TRR(i,j)-1;
 +
            end
 +
           
 +
            CS(i,j+1)= CS(i,j);
 +
            CSa(i,j+1)= CSa(i,j);
 +
            U(i,j+1)= U(i,j);
 +
            O(i,j+1)= O(i,j);
 +
            Of(i,j+1)= Of(i,j);
 +
            Uf(i,j+1)= Uf(i,j);
 +
            TR(i,j+1)= TR(i,j);
 +
            TA(i,j+1)= TA(i,j);
 +
            R(i, j+1)= R(i,j);
 +
           
 +
        end   
 +
    end
 +
end
 +
%
 +
Pro_CS=zeros(numEvents);
 +
Pro_CSa=zeros(numEvents);
 +
Pro_Uf=zeros(numEvents);
 +
Pro_U=zeros(numEvents);
 +
Pro_Of=zeros(numEvents);
 +
Pro_O=zeros(numEvents);
 +
Pro_TR=zeros(numEvents);
 +
Pro_TRR=zeros(numEvents);
 +
Pro_TA=zeros(numEvents);
 +
Pro_R=zeros(numEvents);
 +
%
 +
%
 +
for i=1:numEvents
 +
  Pro_CS(i)=(sum(CS(:,i)))/numCells;
 +
end
 +
for i=1:numEvents
 +
    Pro_CSa(i)=(sum(CSa(:,i)))/numCells;
 +
end
 +
for i=1:numEvents
 +
    Pro_Uf(i)=(sum(Uf(:,i)))/numCells;
 +
end
 +
for i=1:numEvents
 +
    Pro_U(i)=(sum(U(:,i)))/numCells;
 +
end
 +
for i=1:numEvents
 +
  Pro_Of(i)=(sum(Of(:,i)))/numCells;
 +
end
 +
for i=1:numEvents
 +
  Pro_O(i)=(sum(O(:,i)))/numCells;
 +
end
 +
for i=1:numEvents
 +
  Pro_TR(i)=(sum(TR(:,i)))/numCells;
 +
end
 +
for i=1:numEvents
 +
  Pro_TRR(i)=(sum(TRR(:,i)))/numCells;
 +
end
 +
for i=1:numEvents
 +
  Pro_TA(i)=(sum(TA(:,i)))/numCells;
 +
end
 +
for i=1:numEvents
 +
  Pro_R(i)=(sum(R(:,i)))/numCells;
 +
end
 +
%
 +
%
 +
%
 +
function [newTime newData] = regIntervalFixed2(Time,Data,step,timeMax)
 +
%
 +
% This function takes the output of a gillespie's simulation and turns it
 +
% into a matrix with a time vector of regular intervals.
 +
%
 +
% INPUT
 +
% Time: A  matrix sized number of Simulations(Cells) by number of Events
 +
%  that contains the time asociated with each event in each simulation.
 +
% Data: A matrix with the data. (events in columns and simulation in rows)
 +
% step: a scalar of the time interval. It must have the same units as the
 +
%  entrances of the matrix Time.
 +
% timeMax: the maximum time required. It must have the same units as the
 +
%  entrances of the matrix Time.
 +
%
 +
% OUTPUT:
 +
% newTime: a vector of time regularly spaced.
 +
% newData: a matrix where each row represents a simulation (cell) and each
 +
%  column represents the number of molecules for each time in newTime
 +
 
 +
L=size(Time,2);
 +
 
 +
newTime=0:step:timeMax;
 +
newData = zeros(size(Data,1),length(newTime));
 +
 
 +
for k=1:size(Data,1)
 +
    newdata=zeros(size(newTime)); 
 +
    j=1;
 +
    for i=1:length(newTime)
 +
 
 +
      if step*(i-1)>Time(k,j) && j<L
 +
          while step*(i-1)>Time(k,j) && j<L
 +
              j=j+1;
 +
          end   
 +
      end
 +
      newdata(i)=Data(k,j);
 +
    end
 +
    newData(k,:) = newdata;
 +
end
 +
%newdata(L+1)=data(L);
 +
%
 +
%
 +
function impulse= funcImpulso(t)
 +
 
 +
if (t>150 && t<310)
 +
    impulse=100; % Molecules/L
 +
else
 +
    impulse=0; %Molecules/L
 +
end
 +
%
 +
end
 +
%
 +
%
 +
function w=condiciones()
 +
%
 +
xi=[0.5 1000 5 1000];
 +
%
 +
y=fsolve(@CondInd,xi,optimset('algorithm','levenberg-marquardt','maxiter',100000,'tolfun',1e-9));
 +
%
 +
w=y;
 +
%
 +
end
 +
%
 +
%
 +
</textarea>
 +
</p>
 +
</div>
 +
<br>
 +
<br>
 +
<br>
 +
<br>
 +
<br>
 +
 
 +
 
 +
<center>
 +
<div class="button-fill orange" ><div class="button-text"> Back to modeling</div><div class="button-inside"><div class="inside-text"><a style="text-decoration: none; background-color: none; color: red;" href="https://2014.igem.org/Team:Colombia/Modeling">Go! </a></div></div></div>
 +
</center>
 +
<script>
 +
    $(".button-fill").hover(function () {
 +
    $(this).children(".button-inside").addClass('full');
 +
}, function() {
 +
  $(this).children(".button-inside").removeClass('full');
 +
});
 +
    //@ sourceURL=pen.js
 +
  </script>
 +
 
 +
</html>

Latest revision as of 03:10, 18 October 2014

Wheeltz - CSS3 Navigational Wheel Menu

  • Home
  • iGEM
  • Facebook
  • Twitter



Scripting

Feel free to expand and scroll through the text boxes in order to further examine the code.


Deterministic Model
This code creates the differential equations governing the concentration dinamics of each protein in our model, finds the steady state solutions and then solves them using the numerical aproximation method Runge-Kutta


Metropolis-Hastings Algorithm
This code determines which set of values for missing parameters yield the most desirable response.


Sensitivity Analysis
This code points out which are the critical parameters in the system (those that change the response drastically).


Stochastic Model
Sometimes probabilistic models better describe certain systems, This required a fairly big amount of computational power






Back to modeling
Go!

Retrieved from "http://2014.igem.org/Team:Colombia/Scripting"