Team:Colombia/Scripting

From 2014.igem.org

(Difference between revisions)
 
(6 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>
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 readonly>
+
<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 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 readonly>
+
<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 595: Line 594:
<b><font color="#8A0808" size="3" >Stochastic Model</font></b>
<b><font color="#8A0808" size="3" >Stochastic Model</font></b>
<br>
<br>
-
Sometimes probabilistic models better describe certain systems<br>
+
Sometimes probabilistic models better describe certain systems, This required a fairly big amount of computational power<br>
-
<textarea readonly>
+
<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