Team:Colombia/Scripting

From 2014.igem.org

(Difference between revisions)
 
(43 intermediate revisions not shown)
Line 1: Line 1:
-
{{Http://2014.igem.org/Team:Colombia Uniandes/Bootstrap}}
+
{{Http://2014.igem.org/Team:Colombia}}
<html>
<html>
<div class="span11" style="text-align: justify;">
<div class="span11" style="text-align: justify;">
-
<br><br>
+
<br>
<p>
<p>
-
I smiled and shook my head. "I can quite understand your thinking so." I said. "Of course, in your position of unofficial adviser and helper to everybody who is absolutely puzzled, throughout three continents, you are brought in contact with all that is strange and bizarre. But here"--I picked up the morning paper from the ground--"let us put it to a practical test. Here is the first heading upon which I come. 'A husband's cruelty to his wife.' There is half a column of print, but I know without reading it that it is all perfectly familiar to me. There is, of course, the other woman, the drink, the push, the blow, the bruise, the sympathetic sister or landlady. The crudest of writers could invent nothing more crude."
+
<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
 +
%
-
"Indeed, your example is an unfortunate one for your argument," said Holmes, taking the paper and glancing his eye down it. "This is the Dundas separation case, and, as it happens, I was engaged in clearing up some small points in connection with it. The husband was a teetotaler, there was no other woman, and the conduct complained of was that he had drifted into the habit of winding up every meal by taking out his false teeth and hurling them at his wife, which, you will allow, is not an action likely to occur to the imagination of the average story-teller. Take a pinch of snuff, Doctor, and acknowledge that I have scored over you in your example."
+
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
-
He held out his snuffbox of old gold, with a great amethyst in the centre of the lid. Its splendour was in such contrast to his homely ways and simple life that I could not help commenting upon it.
+
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);                                                            %
-
"Ah," said he, "I forgot that I had not seen you for some weeks. It is a little souvenir from the King of Bohemia in return for my assistance in the case of the Irene Adler papers."
+
%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
-
"And the ring?" I asked, glancing at a remarkable brilliant which sparkled upon his finger.
+
L=size(Time,2);
-
"It was from the reigning family of Holland, though the matter in which I served them was of such delicacy that I cannot confide it even to you, who have been good enough to chronicle one or two of my little problems."
+
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"