Team:Colombia/Scripting

From 2014.igem.org

(Difference between revisions)
 
(40 intermediate revisions not shown)
Line 2: Line 2:
<html>
<html>
<div class="span11" style="text-align: justify;">
<div class="span11" style="text-align: justify;">
-
<br><br>
+
<br>
<p>
<p>
-
<script>
+
<br>
-
function y=CondIni(x)
+
<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
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                 %
+
%                               VARIABLES                                 %
%--------------------------------------------%
%--------------------------------------------%
%
%
-
  C=0;                                           % Extracellular concentration of the cholerae autoinducer-1 (CAI-1)
+
  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)
  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)
+
  CSa=x(2);                                     % Concentration of the membrane bound CqsS protein (CAI-1 bound=Active)
  Uf= x(3);                                      % Phosphorelay protein LuxU (phosphorylated) Concentration
  Uf= x(3);                                      % Phosphorelay protein LuxU (phosphorylated) Concentration
  U=x(4);                                        % Phosphorelay protein LuxU (unphosphorylated) Concentration
  U=x(4);                                        % Phosphorelay protein LuxU (unphosphorylated) Concentration
Line 23: Line 107:
  TR=x(7);                                      % Ptet Repressor protein concentration
  TR=x(7);                                      % Ptet Repressor protein concentration
  TA=x(8);                                      % Ptet Activator protein concentration
  TA=x(8);                                      % Ptet Activator protein concentration
-
  R=x(9);                                       % Response molecule concentration
+
  R=x(9);                                         % Response molecule concentration
%  
%  
%--------------------------------------------%
%--------------------------------------------%
-
%                 PARAMETERS                 %
+
%                               PARAMETERS                             %
%--------------------------------------------%
%--------------------------------------------%
%
%
Line 32: Line 116:
kcd=1;                                          % CAI1 and CqsS decoupling Rate
kcd=1;                                          % CAI1 and CqsS decoupling Rate
dsu=4;                                          % LuxU* dephosphorylation rate through CqsS*
dsu=4;                                          % LuxU* dephosphorylation rate through CqsS*
-
duo=4;                                         % LuxO* dephosphorylation rate through LuxU
+
duo=4;                                         % LuxO* dephosphorylation rate through LuxU
fsu=2;                                          % LuxU* phosphorylation rate through CqsS
fsu=2;                                          % LuxU* phosphorylation rate through CqsS
fuo=2;                                          % LuxO* phosphorylation rate through LuxU*
fuo=2;                                          % LuxO* phosphorylation rate through LuxU*
%
%
gcs=1;                                          % CqsS protein decay rate
gcs=1;                                          % CqsS protein decay rate
-
gcsa=1;                                         % CqsS* protein decay rate
+
gcsa=1;                                       % CqsS* protein decay rate
guf=1;                                          % LuxU* protein decay rate
guf=1;                                          % LuxU* protein decay rate
gu=1;                                          % LuxU protein decay rate
gu=1;                                          % LuxU protein decay rate
Line 49: Line 133:
au=3;                                          % LuxU basal production rate
au=3;                                          % LuxU basal production rate
ao=3;                                          % LuxO basal production rate
ao=3;                                          % LuxO basal production rate
-
ar=0.01;                                       % response molecule basal production rate
+
ar=0.01;                                       % response molecule basal production rate
-
atr=0.01;                                       % TR basal production rate
+
atr=0.01;                                     % TR basal production rate
-
ata=0.01;                                       % TA basal production rate
+
ata=0.01;                                     % TA basal production rate
-
%
+
 
btr=5;                                          % Maximum rate of TR expression
btr=5;                                          % Maximum rate of TR expression
bta=5;                                          % Maximum rate of TA expression
bta=5;                                          % Maximum rate of TA expression
br=5;                                          % Maximum rate of response molecule expression
br=5;                                          % Maximum rate of response molecule expression
-
ho=1.5;                                         % LuxO*- DNA coupling rate
+
ho=1.5;                                       % LuxO*- DNA coupling rate
htr=2;                                          % TRdomain-DNA coupling rate
htr=2;                                          % TRdomain-DNA coupling rate
%
%
Line 62: Line 146:
%
%
%--------------------------------------------%
%--------------------------------------------%
-
%                 Equations                 %
+
%                                 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
+
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
+
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
+
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
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
+
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
+
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
+
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
+
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
+
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(1)=dCS;                                      %
-
y(2)=dCSa;                                     %
+
y(2)=dCSa;                                     %
y(3)=dUf;                                      %
y(3)=dUf;                                      %
y(4)=dU;                                        %
y(4)=dU;                                        %
Line 83: Line 167:
y(7)=dTR;                                      %
y(7)=dTR;                                      %
y(8)=dTA;                                      %
y(8)=dTA;                                      %
-
y(9)=dR;                                       %
+
y(9)=dR;                                         %
%
%
-
y=y';                                           % Transposing the vector because of matlab language restrictions
+
y=y';                                             % Transposing the vector because of matlab language restrictions
%
%
end
end
%
%
-
</script>
+
%
 +
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>
</p>
</div>
</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>
</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"