Team:Colombia/Scripting

From 2014.igem.org

(Difference between revisions)
 
(14 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>
-
<br>
 
<br>
<br>
<center><b> <font size="24">
<center><b> <font size="24">
-
<br><br>
+
<br>
Scripting
Scripting
</font><br>
</font><br>
-
<h1  class="curs1"><font size="6">
+
<h1  class="curs1"><font size="4">
Feel free to expand and scroll through the text boxes in order to further examine the code.
Feel free to expand and scroll through the text boxes in order to further examine the code.
</font></h1></b></center>
</font></h1></b></center>
Line 17: Line 16:
<br>
<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>
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>
+
<textarea readonly rows = "20">
function y=det()
function y=det()
%
%
Line 194: Line 193:
<br>
<br>
This code determines which set of values for missing parameters yield the most desirable response.<br>
This code determines which set of values for missing parameters yield the most desirable response.<br>
-
<textarea>
+
<textarea readonly rows = "20">
function y=model(patr, pbtr, pho, pkttr, pata, pbta, phtr,pn)
function y=model(patr, pbtr, pho, pkttr, pata, pbta, phtr,pn)
%
%
Line 201: Line 200:
%  
%  
%--------------------------------------------%
%--------------------------------------------%
-
%                 PARAMETERS                 %
+
%                             PARAMETERS                             %
%--------------------------------------------%
%--------------------------------------------%
%
%
Line 247: Line 246:
%
%
%--------------------------------------------%
%--------------------------------------------%
-
%     Runge Kutta 4 order aproximation     %
+
%               Runge Kutta 4 order aproximation           %
%--------------------------------------------%
%--------------------------------------------%
%
%
Line 304: Line 303:
%xlabel('Time [h]');                                                %
%xlabel('Time [h]');                                                %
%ylabel('Concentration [nanomolar]');                    %
%ylabel('Concentration [nanomolar]');                    %
-
%title('Response');                                                 %
+
%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)];
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)];
Line 441: Line 440:
function y=input_Metropolis(t)
function y=input_Metropolis(t)
%--------------------------------------------%
%--------------------------------------------%
-
%         CAI-1 (C) pulse simulation         %
+
%                   CAI-1 (C) pulse simulation                   %
%--------------------------------------------%
%--------------------------------------------%
if t>100 && t<200   
if t>100 && t<200   
Line 454: Line 453:
<br>
<br>
This code points out which are the critical parameters in the system (those that change the response drastically).<br>
This code points out which are the critical parameters in the system (those that change the response drastically).<br>
-
<textarea>
+
<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
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
Line 460: Line 459:
%  
%  
%--------------------------------------------%
%--------------------------------------------%
-
%                 PARAMETERS                 %
+
%                               PARAMETERS                           %
%--------------------------------------------%
%--------------------------------------------%
%
%
-
kcc=1;                                             % CAI1 and CqsS coupling Rate
+
kcc=1;                                                           % CAI1 and CqsS coupling Rate
-
kcd=1;                                             % CAI1 and CqsS decoupling Rate
+
kcd=1;                                                           % CAI1 and CqsS decoupling Rate
-
kcu=2;                                             % Michaelis Menten Constant for the phoshporylation CqsS-LuxU
+
kcu=2;                                                           % Michaelis Menten Constant for the phoshporylation CqsS-LuxU
-
kuc=1;                                             % Michaelis Menten Constant of the phosphase process CqSS-LuxU
+
kuc=1;                                                           % Michaelis Menten Constant of the phosphase process CqSS-LuxU
-
kuo=2;                                             % Michaelis Menten Constant for the phoshporylation LuxU-LuxO
+
kuo=2;                                                           % Michaelis Menten Constant for the phoshporylation LuxU-LuxO
-
kou=1;                                             % Michaelis Menten Constant of the phosphase process LuxU-LuxO
+
kou=1;                                                           % Michaelis Menten Constant of the phosphase process LuxU-LuxO
-
kttr=5;                                             % Dimerization rate of TetR ---[Sensitivity Analysis]--- [5-50/5]
+
kttr=5;                                                           % Dimerization rate of TetR ---[Sensitivity Analysis]--- [5-50/5]
%
%
-
gcs=0.5;                                           % CqsS protein decay rate
+
gcs=0.5;                                                         % CqsS protein decay rate
-
gcsa=1;                                             % CqsS* protein decay rate
+
gcsa=1;                                                         % CqsS* protein decay rate
-
guf=0.12;                                           % LuxU* protein decay rate
+
guf=0.12;                                                       % LuxU* protein decay rate
-
gu=0.65;                                           % LuxU protein decay rate
+
gu=0.65;                                                       % LuxU protein decay rate
-
gof=0.12;                                           % LuxO* protein decay rate
+
gof=0.12;                                                       % LuxO* protein decay rate
-
go=0.12;                                           % LuxO protein decay rate
+
go=0.12;                                                       % LuxO protein decay rate
-
gtr=0.12;                                           % Ptet Repressor protein decay rate
+
gtr=0.12;                                                       % Ptet Repressor protein decay rate
-
gttr=0.12;                                         % Ptet2 dimer Repressor decay rate
+
gttr=0.12;                                                     % Ptet2 dimer Repressor decay rate
-
gta=0.12;                                           % Ptet Activator protein decay rate
+
gta=0.12;                                                       % Ptet Activator protein decay rate
-
gr=0.12;                                           % Response molecule decay rate
+
gr=0.12;                                                         % Response molecule decay rate
%
%
-
acs=5;                                             % CS basal production rate
+
acs=5;                                                             % CS basal production rate
-
au=5;                                               % LuxU basal production rate
+
au=5;                                                             % LuxU basal production rate
-
ao=5;                                               % LuxO basal production rate
+
ao=5;                                                             % LuxO basal production rate
-
ar=0;                                               % response molecule basal production rate
+
ar=0;                                                               % response molecule basal production rate
-
atr=1;                                             % TR basal production rate  ---[Sensitivity Analysis]--- [0.5-5/0.5]
+
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]
+
ata=0.1;                                                           % TA basal production rate ---[Sensitivity Analysis]--- [0.1-10/0.5]
%
%
-
bcu=1;                                             % Phosphorylation rate of CqsS to LuxU
+
bcu=1;                                                             % Phosphorylation rate of CqsS to LuxU
-
buc=0.1;                                           % Phosphase rate of CqsS to LuxU
+
buc=0.1;                                                         % Phosphase rate of CqsS to LuxU
-
buo=3.2;                                           % Phosphorylation rate LuxU-LuxO
+
buo=3.2;                                                         % Phosphorylation rate LuxU-LuxO
-
bou=1.58;                                           % Phosphase rate LuxU-LuxO
+
bou=1.58;                                                       % Phosphase rate LuxU-LuxO
-
btr=5;                                             % Maximum rate of TR expression ---[Sensitivity Analysis]--- [1-10/1]
+
btr=5;                                                             % Maximum rate of TR expression ---[Sensitivity Analysis]--- [1-10/1]
-
bta=5;                                             % Maximum rate of TA 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]
+
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]
+
htr=5;                                                               % TRdomain-DNA coupling rate ---[Sensitivity Analysis]--- [1-10/1]
%
%
-
n=2;                                               % Hill's coefficient ---[Sensitivity Analysis]--- [1-4/0.5]
+
n=2;                                                                 % Hill's coefficient ---[Sensitivity Analysis]--- [1-4/0.5]
%
%
-
h=400;                                             % Maximum time
+
h=400;                                                             % Maximum time
-
m=1/10;                                             % Step lenght [h]
+
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.
+
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.
+
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
+
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)
+
     kttr= param(w);                                               % New parameter value (This MUST be changed everytime a parameter is being evaluated)
%
%
%
%
%--------------------------------------------%
%--------------------------------------------%
-
%     Runge Kutta 4 order aproximation     %
+
%               Runge Kutta 4 order aproximation             %
%--------------------------------------------%
%--------------------------------------------%
%
%
%
%
%
%
-
conInd=[1,0,1,1,1,1,0,1,1,1];                       %
+
conInd=[1,0,1,1,1,1,0,1,1,1];                               %
-
l=(0:m:h)';                                         % Time vector
+
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
+
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
+
C=zeros(1,length(l));                                             % zero vector to be filled
-
x(1,:)=conInd;                                     %
+
x(1,:)=conInd;                                                       %
%
%
for k=1:length(l)-1
for k=1:length(l)-1
-
     xk=x(k,:);                                     % Capture of the last value in the matrix, the actual variable values
+
     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
+
     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
+
     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
+
     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
+
     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/6*(k1+2*k2+2*k3+k4)';                  % New value calculation for the variables
-
     xk2=zeros(1,length(xk1));                       %
+
     xk2=zeros(1,length(xk1));                                 %
%     
%     
     for p=1:length(xk1)
     for p=1:length(xk1)
Line 539: Line 538:
         end   
         end   
     end
     end
-
     x(k+1,:)=xk2;                       % Actualizacion del nuevo vector de variables en la matriz
+
     x(k+1,:)=xk2;                               % Actualizacion del nuevo vector de variables en la matriz
end
end
%
%
Line 551: Line 550:
%
%
%
%
-
   if (R(it)-R(1011))/(R(1011))> 1.5     % If the response changes 1.5-fold relative to previous iteration
+
   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;
   rta(1,w)= l(it)-100;
Line 590: Line 589:
ylabel('Change in response (nM)')
ylabel('Change in response (nM)')
title('Sensitivity analysis')
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>
</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