Team:Colombia Uniandes/Scripting

From 2013.igem.org

(Difference between revisions)
(→Stochastic)
(→Equation solver)
Line 572: Line 572:
====Equation solver====
====Equation solver====
-
 
+
%
-
 
+
%
-
clear all
+
clear all
-
clc
+
clc
-
 
+
%
-
 
+
%
-
 
+
%
  %---------Parameters------%
  %---------Parameters------%
-
global gammaN Kp beta Kd Kx alfaR deltaR Kt deltaT n alfaP deltaP
+
%
-
 
+
%
-
gammaN=0.5034e-4; %Diffussion rate of Nickel (1/min)
+
%
-
Kp=0.000634; %Dynamic constant for the entrance of nickel to the cell
+
global gammaN Kp beta Kd Kx alfaR deltaR Kt deltaT n alfaP deltaP
-
beta=0.166; %Porine maximum expression rate (nM/min)
+
%
-
Kd=276e-3; %Association constant DNA and repressor (nM)
+
gammaN=0.5034e-4; %Diffussion rate of Nickel (1/min)
-
Kx=25e-3; %Association constant of the repressor with nickel (nM)
+
%
-
alfaR= 5; %Represor basal production rate (nM/min)
+
Kp=0.000634; %Dynamic constant for the entrance of nickel to the cell
-
deltaR=1/1200; % Represor destruction rate (1/min)
+
%
-
Kt=820e-3;% Rate constant for the formation of the tetramer (Giraldo et al)
+
beta=0.166; %Porine maximum expression rate (nM/min)
-
deltaT=1/1200; %Tetramer destruction rate (!/min)  
+
%
-
n=1; %Hill coefficient (cooperation constant)
+
Kd=276e-3; %Association constant DNA and repressor (nM)
-
alfaP=0.031; %Porine basal production rate (nM/min)
+
%
-
deltaP=1/1200;%Porine destruction rate (1/min)
+
Kx=25e-3; %Association constant of the repressor with nickel (nM)
-
   
+
%
-
 
+
alfaR= 5; %Represor basal production rate (nM/min)
-
yo=[0 0 0 0];
+
%
-
con=fsolve(@(x)CondIniciales(x),yo, optimset('display','iter','MaxIter',1000000,'algorithm','levenberg-marquardt','tolfun',1e-9));
+
deltaR=1/1200; % Represor destruction rate (1/min)
-
assignin('base','ini',con);
+
%
-
 
+
Kt=820e-3;% Rate constant for the formation of the tetramer (Giraldo et al)
-
cond=abs(con);
+
%
-
 
+
deltaT=1/1200; %Tetramer destruction rate (!/min)
-
h=30; %Tiempo maximo
+
%
-
m=0.01; %Paso
+
n=1; %Hill coefficient (cooperation constant)
-
l=(0:m:h);%Vector de tiemp
+
%
-
 
+
alfaP=0.031; %Porine basal production rate (nM/min)
-
condI=[9.88e3 cond];
+
%
-
x=zeros(length(l),length(condI));
+
deltaP=1/1200;%Porine destruction rate (1/min)
-
x(1,:)=condI;
+
  %
-
 
+
%
-
for k=1:length(l)-1
+
%
-
   
+
yo=[0 0 0 0];
-
     xk=x(k,:); %Captura de la ultima posicion de la matirz, es decir, los
+
%
-
    %valores actuales de las variables
+
con=fsolve(@(x)CondIniciales(x),yo, optimset('display','iter','MaxIter',1000000,'algorithm','levenberg-marquardt','tolfun',1e-9));
-
   
+
%
 +
assignin('base','ini',con);
 +
%
 +
%
 +
cond=abs(con);
 +
%
 +
h=30; %Tiempo maximo
 +
%
 +
m=0.01; %Paso
 +
l=(0:m:h);%Vector de tiemp
 +
%
 +
condI=[9.88e3 cond];
 +
x=zeros(length(l),length(condI));
 +
x(1,:)=condI;
 +
%
 +
for k=1:length(l)-1
 +
%   
 +
     xk=x(k,:); %Captura de la ultima posicion de la matirz, es decir, los valores actuales de las variables
 +
% 
     k1=EqNick(xk,l(k)); %Primera pendiente del metodo de RK4
     k1=EqNick(xk,l(k)); %Primera pendiente del metodo de RK4
     k2=EqNick(xk+(m/2*k1)',l(k)+m/2); %Segunda pendiente del metodo de RK4
     k2=EqNick(xk+(m/2*k1)',l(k)+m/2); %Segunda pendiente del metodo de RK4
     k3=EqNick(xk+(m/2*k2)',l(k)+m/2); %Tercera pendiente del metodo de RK4
     k3=EqNick(xk+(m/2*k2)',l(k)+m/2); %Tercera pendiente del metodo de RK4
     k4=EqNick(xk+(m*k3)',l(k)+m); %Cuarta pendiente del metodo de RK4
     k4=EqNick(xk+(m*k3)',l(k)+m); %Cuarta pendiente del metodo de RK4
-
   
+
% 
     xk1=xk+m/6*(k1+2*k2+2*k3+k4)'; %Calculo de nuevos valores para las
     xk1=xk+m/6*(k1+2*k2+2*k3+k4)'; %Calculo de nuevos valores para las
     %variables
     %variables
-
   
+
% 
     %xk1=xk+m*ecuaDif(l(k),xk)'; %Method of Newton
     %xk1=xk+m*ecuaDif(l(k),xk)'; %Method of Newton
-
   
+
% 
     xk2=zeros(1,length(xk1));
     xk2=zeros(1,length(xk1));
-
   
+
% 
-
   
+
% 
     for p=1:length(xk1)
     for p=1:length(xk1)
-
       
+
%     
         if(xk1(p)<0.00000001)
         if(xk1(p)<0.00000001)
-
           
+
%         
             xk2(p)=0;
             xk2(p)=0;
         else
         else
-
           
+
%         
             xk2(p)=xk1(p);
             xk2(p)=xk1(p);
         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
-
 
+
%
-
No=x(:,1);
+
No=x(:,1);
-
Ni=x(:,2);
+
Ni=x(:,2);
-
assignin('base','No',No);
+
assignin('base','No',No);
-
assignin('base','Nistable',Ni(length(Ni)));
+
assignin('base','Nistable',Ni(length(Ni)));
-
disp(Ni(length(Ni)))
+
disp(Ni(length(Ni)))
-
R=x(:,3);
+
R=x(:,3);
-
assignin('base','Rstable',R(length(R)));
+
assignin('base','Rstable',R(length(R)));
-
T=x(:,4);
+
T=x(:,4);
-
assignin('base','Tstable',T(length(T)));
+
assignin('base','Tstable',T(length(T)));
-
P=x(:,5);
+
P=x(:,5);
-
assignin('base','Pstable',P(length(P)));
+
assignin('base','Pstable',P(length(P)));
-
cond=[R(length(R)) T(length(T)) P(length(P))];
+
cond=[R(length(R)) T(length(T)) P(length(P))];
-
assignin('base','cond',cond);
+
assignin('base','cond',cond);
-
figure(1)
+
figure(1)
-
plot(l,No,l,P)
+
plot(l,No,l,P)
-
legend('No','P')
+
legend('No','P')
-
xlabel('Time (min)')
+
xlabel('Time (min)')
-
ylabel('Concentration (nM)')
+
ylabel('Concentration (nM)')
-
figure(2)
+
figure(2)
-
plot(l,Ni)
+
plot(l,Ni)
-
legend('Ni')
+
legend('Ni')
-
xlabel('Time (min)')
+
xlabel('Time (min)')
-
ylabel('Concentration (nM)')
+
ylabel('Concentration (nM)')
-
figure(3)
+
figure(3)
-
plot(l,R)
+
plot(l,R)
-
legend('R')
+
legend('R')
-
xlabel('Time (min)')
+
xlabel('Time (min)')
-
ylabel('Concentration (nM)')
+
ylabel('Concentration (nM)')
-
figure(4)
+
figure(4)
-
plot(l,T)
+
plot(l,T)
-
legend('T')
+
legend('T')
-
xlabel('Time (min)')
+
xlabel('Time (min)')
-
ylabel('Concentration (nM)')
+
ylabel('Concentration (nM)')
-
figure(5)
+
figure(5)
-
plot(l,P)
+
plot(l,P)
-
legend('P')
+
legend('P')
-
xlabel('Time (min)')
+
xlabel('Time (min)')
-
ylabel('Concentration (nM)')
+
ylabel('Concentration (nM)')  
 +
%
===Stochastic===
===Stochastic===

Revision as of 22:30, 25 September 2013


Scripting

Contents

Glucocorticoid Detection System


Deterministic model

Equations

function y = EcuacionesGluco(t,x)
global gammaGR mGRIR mCC deltaGRI alfaR deltaR deltaCC betaCC k n deltaS H
%
%---------Parameters------%
%
%
GRO=funcImpulso(t);
%
%------ Variables%------%
%
GRI= x(1); %Glucocorticoid inside the cell 
R=x(2); %Receptor in the cytoplasm
CC=x(3); %Receptor -Glucocorticoid complex
V=x(4); %Violacein
%
%
%---Equations---%
dGRI=gammaGR*(GRO-GRI)-mGRIR*GRI*R+mCC*CC-deltaGRI*GRI;
dR=alfaR-mGRIR*GRI*R+mCC*CC-deltaR*R;
dCC=mGRIR*GRI*R-mCC*CC-deltaCC*CC-(betaCC*CC.^n)/(k^n+CC.^n);%Revisar
dV=H*(betaCC*CC^n)/(k^n+CC^n)-deltaS*V;
%
y1(1)=dGRI;
y1(2)=dR;
y1(3)=dCC;
y1(4)=dV;
%
%
y= y1';
%
end

Equation solver

%global gammaGR mGRIR mCC deltaGRI alfaR deltaR deltaCC betaCC k n deltaS H
%
%
gammaGR= 0.1;    %Diffussion rate of glucocorticoid inside the cell (mm/min)
%
mGRIR=1.080e-3;        % GRI-R complex formation kinetic constant (1/umol min)
%
mCC=1.14*10^-8;    %GRI-R Complex formation reverse  kinetic constant (1/min)
%
deltaGRI=0.00833;  %Glucocorticoids Destruction rate inside the cell (1/min)
%
alfaR= 0.8e3;           %Basal production rate of the receptor (umol/min) 
%
deltaR=0.004166;          %Receptor destruction rate inside the cell (1/min)
deltaCC=0.004166;  % GRI-R complex Destruction rate (1/min)
betaCC=0.5e3;     % GRI-R complex maximum expression rate (umol/min)
k=0.05e3;          %Hill's constant for the GRI-R complex dimmer binding to his respective region (umol)
n=2;               %Hill coefficient (cooperation constant)
deltaS=0.04166;   %Signal destruction rate (1/min)
H=2;               %Correction constant for the signal
%
%
%
h=60; %Tiempo maximo
%
m=0.01; %Longitud de paso [s]
%
t=0:m:h; %Vector tiempo
%
xi=[0 0 0 0];
%
y=fsolve(@CondIndGluco,xi,optimset('algorithm','levenberg-marquardt','maxiter',100000,'tolfun',1e-9));
%
conInd=y;
assignin('base','conInd',conInd);
l=(0:m:h)'; %Vector de tiempo
%
x=zeros(length(l),length(conInd)); %Matriz de variables, en las columnas varia
%la variable y en las filas varia el tiempo
%
GRO=zeros(1,length(l));
%
x(1,:)=conInd;
%
for u=1:length(l)-1
%    
   xk=x(u,:); %Captura de la ultima posicion de la matirz, es decir, los
   %valores actuales de las variables
%    
   k1=EcuacionesGluco(l(u),xk); %Primera pendiente del metodo de RK4
   k2=EcuacionesGluco(l(u)+m/2,xk+(m/2*k1)'); %Segunda pendiente del metodo de RK4
   k3=EcuacionesGluco(l(u)+m/2,xk+(m/2*k2)'); %Tercera pendiente del metodo de RK4
   k4=EcuacionesGluco(l(u)+m,xk+(m*k3)'); %Cuarta pendiente del metodo de RK4
%   
   xk1=xk+m/6*(k1+2*k2+2*k3+k4)'; %Calculo de nuevos valores para las
   %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(u+1,:)=xk2; %Actualizacion del nuevo vector de variables en la matriz
%   
%   
%   
%   
%   
end
%
for j=1:length(l)
%   
   if (l(j)<(10) || l(j)>(30))
%       
       GRO(j)=155;
%       
   else
%       
       GRO(j)=155*1.3;
%       
%       
   end
%   
%   
end
%
GRI=x(:,1);
R=x(:,2);
CC=x(:,3);
V=x(:,4);
%
%
figure(1) 
plot(l,R)%,l,GRO)%,l,CC,l,V)
legend('Receptor')%,'Glucocorticoid') %, 'Complex', 'Signal')
xlabel('Time')
ylabel('Concetration (micromolar)')
title('Glucocorticoid model')
%
figure(2)
plot(l,CC)%,l,GRO)
legend('Complejo')%,'Glucocorticoid')
%
figure(3)
plot(l,V)%,l,GRO)
legend('Senal')%,'Glucocorticoid')
%
figure(4)
plot(l,GRI)%,l,GRO)
legend('GRI')%,'Glucocorticoid')
%

Stochastic

%Stochastic Simulation of the Chemical Reactions of the Glucocorticoid 

   %System
   %

   %%----Parameters of the System----
   global gammaGR mGRIR mCC deltaGRI alfaR deltaR deltaCC betaCC k n deltaS H

   %%-----Initial Conditions of the glucocorticoid system-----
   %Run CondIndGluco.m with the parameters of the system.

   con=condiciones();

   GRO_0= round(10000); %Glucocorticoid Outside the Cell
   GRI_0 = round(con(1)); %Glucocorticoid Inside the Cell
   R_0= round(con(2)); %Receptor
   CC_0= round(con(3)); %Receptor- Glucocorticoid Complex
   V_0= round(con(4));%Signal

   %%----Simulation Constants----
   numCells = 10;
   numEvents= 60000;
   numCurves = 5; % Restriction(numCurves<numCells)
   Step=0.1;
   tMax=45; %Minutes

   %Plot the behaviour of the molecule 1/0 for yes/no
   Plot_GRO = 1;            
   Plot_GRI=1;            
   Plot_R =1;             
   Plot_CC = 1;             
   Plot_V = 1;   
   PlotMean=1;          

   %%----Stochastic Simulation----
   %Variables of the System
   GRO= zeros(numCells,numEvents);
   GRI= zeros(numCells,numEvents);
   R= zeros(numCells,numEvents);
   CC= zeros(numCells,numEvents);
   V = zeros(numCells,numEvents);

   t = zeros(numCells,numEvents); %Time
   tr = zeros(numCells,numEvents);
   %Initial Conditions
   GRO(:,1)=GRO_0*ones(numCells,1);
   GRI(:,1)=GRI_0*ones(numCells,1);
   R(:,1)=R_0*ones(numCells,1);
   CC(:,1)=CC_0*ones(numCells,1);
   V(:,1)=V_0*ones(numCells,1);

   %Stochastic Simulation

   for S_ns =1:numCells



       for S_ne = 1:numEvents
           %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= gammaGR*GRO(S_ns,S_ne);
           e2= gammaGR*GRI(S_ns,S_ne);
           e3= mGRIR*GRI(S_ns,S_ne)*R(S_ns,S_ne);
           e4= mCC*CC(S_ns,S_ne);
           e5= deltaGRI*GRI(S_ns,S_ne);
           e6= alfaR;
           e7= deltaR*R(S_ns,S_ne);
           e8= deltaCC*CC(S_ns,S_ne);
           %e8= (betaCC*((CC(S_ns,S_ne))^n))/(k^n)+ ((CC(S_ns,S_ne))^n);
           e9= H*((betaCC*((CC(S_ns,S_ne))^n))/((k^n)+ ((CC(S_ns,S_ne))^n)));
           e10= deltaS*V(S_ns,S_ne);

           %Vector of Events
           vecEvents = [e1,e2,e3,e4,e5,e6,e7,e8,e9,e10];

           %Sum of Events
           sumEvents = sum(vecEvents);

           %Waitiing time until the next ocurrence of a reaction time
           %Probabilistic Distribution
           time= (1/sumEvents)*log(1/rand1);
           tr(S_ns,S_ne)=time;
           t(S_ns,S_ne+1)= t(S_ns,S_ne)+ time;



           %Glucocorticoid outside the cell        

           GRO(S_ns,S_ne+1)=funcImpulso(t(S_ns,S_ne+1));
           if(t(S_ns,S_ne)<130)
               GRI(S_ns,S_ne)=funcImpulso(t(S_ns,S_ne+1))*957/2575;
           end

           %Simulation Creation/Destruction of each molecule
           if rand2<(vecEvents(1)/sumEvents)

               GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne)+1;
               R(S_ns, S_ne+1)= R(S_ns,S_ne);
               CC(S_ns, S_ne+1)= CC(S_ns,S_ne);
               V(S_ns, S_ne+1) = V(S_ns,S_ne);

           elseif rand2< (sum(vecEvents(1:2)))/sumEvents
               if GRI(S_ns,S_ne)<=0

                   GRI(S_ns,S_ne+1)= 0;
               else
                   GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne)-1;
               end

               R(S_ns, S_ne+1)= R(S_ns,S_ne);
               CC(S_ns, S_ne+1)= CC(S_ns,S_ne);
               V(S_ns, S_ne+1) = V(S_ns,S_ne);



           elseif rand2< (sum(vecEvents(1:3)))/sumEvents

               if R(S_ns,S_ne)<=0 || GRI(S_ns,S_ne)<=0
                  CC(S_ns, S_ne+1)= CC(S_ns,S_ne);
                  R(S_ns, S_ne+1)= R(S_ns,S_ne);
                  GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne);
               else
                  CC(S_ns, S_ne+1)= CC(S_ns,S_ne)+1;
                  R(S_ns, S_ne+1)= R(S_ns,S_ne)-1;
                  GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne)-1;
               end
               V(S_ns, S_ne+1) = V(S_ns,S_ne);



           elseif rand2< (sum(vecEvents(1:4)))/sumEvents

               if CC(S_ns,S_ne)<=50
                  CC(S_ns,S_ne+1)=0;
                  GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne);
                  R(S_ns, S_ne+1)= R(S_ns,S_ne); 
               else
                  CC(S_ns, S_ne+1)= CC(S_ns,S_ne)-50;
                  GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne)+50;
                  R(S_ns, S_ne+1)= R(S_ns,S_ne)+50; 
               end

               V(S_ns, S_ne+1) = V(S_ns,S_ne);

           elseif rand2< (sum(vecEvents(1:5)))/sumEvents

               if GRI(S_ns,S_ne)<=0
                  GRI(S_ns,S_ne+1)=0;
               else
                  GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne)-10;
               end



               R(S_ns, S_ne+1)= R(S_ns,S_ne);
               CC(S_ns, S_ne+1)= CC(S_ns,S_ne);
               V(S_ns, S_ne+1) = V(S_ns,S_ne); 

           elseif rand2< (sum(vecEvents(1:6)))/sumEvents

               GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne);
               R(S_ns, S_ne+1)= R(S_ns,S_ne)+600;
               CC(S_ns, S_ne+1)= CC(S_ns,S_ne);
               V(S_ns, S_ne+1) = V(S_ns,S_ne);

           elseif rand2< (sum(vecEvents(1:7)))/sumEvents

                if R(S_ns,S_ne)<=0
                      R(S_ns,S_ne+1)=0;
               else
                   R(S_ns, S_ne+1)= R(S_ns,S_ne)-300;
               end



               GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne);
               CC(S_ns, S_ne+1)= CC(S_ns,S_ne);
               V(S_ns, S_ne+1) = V(S_ns,S_ne);

           elseif rand2< (sum(vecEvents(1:8)))/sumEvents
               if CC(S_ns,S_ne)<=35
                  CC(S_ns,S_ne+1)=round(CC(S_ns,S_ne)*0.75);
               else
                 CC(S_ns, S_ne+1)= CC(S_ns,S_ne)-35;

               end
               GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne);
               R(S_ns, S_ne+1)= R(S_ns,S_ne);                   
               V(S_ns, S_ne+1) = V(S_ns,S_ne);

           elseif rand2< (sum(vecEvents(1:9)))/sumEvents

               GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne);
               R(S_ns, S_ne+1)= R(S_ns,S_ne);
               CC(S_ns, S_ne+1)= CC(S_ns,S_ne);
               V(S_ns, S_ne+1) = V(S_ns,S_ne)+1000;

           else
               if V(S_ns,S_ne)<=0
                  V(S_ns,S_ne+1)=0;
               else
                  V(S_ns, S_ne+1)= V(S_ns,S_ne)-500;
               end

               GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne);
               R(S_ns, S_ne+1)= R(S_ns,S_ne);
               CC(S_ns, S_ne+1)= CC(S_ns,S_ne);

           end
       end
   end

   %% ------------- Time Regularization -----------------------------
   time=t;

   [orgTime, newGRO] = regIntervalFixed2(time, GRO, Step, tMax);
   [orgTime, newGRI] = regIntervalFixed2(time,GRI, Step, tMax);
   [orgTime, newR] = regIntervalFixed2(time, R, Step, tMax);
   [orgTime, newCC] = regIntervalFixed2(time, CC, Step, tMax);
   [orgTime, newV] = regIntervalFixed2(time, V, Step, tMax);

   GROMean=mean(newGRO);
   GRIMean=mean(newGRI);
   RMean=mean(newR);
   CCMean=mean(newCC);
   VMean=mean(newV);



   %% Plotting 
   % figure(1)
   % plot(t(1,:),GRI(1,:),t(1,:),GRO)
   % title('GRI')
   % figure(2)
   % plot(t(1,:),R(1,:))
   % title('R')
   % figure(3)
   % plot(t(1,:),V(1,:))
   % title('V')
   % figure(4)
   % plot(t(1,:),CC(1,:))
   % title('CC')

   if Plot_GRO
       figure,

       for i= 1:numCurves
           stairs(orgTime,newGRO(i,:))
           hold on
       end
       %if PlotMean
           %stairs(orgTime,GROMean,'r')
       %end
       hold off 

       title();
       xlabel('Time(min)');
       ylabel('Glucocorticoid Outside The Cell(Molecules)');
   end

   if Plot_GRI
       figure,
       for i= 1:numCurves
           stairs(orgTime,newGRI(i,:))

           hold on
       end

       if PlotMean
           stairs(orgTime,GRIMean,'r','LineWidth',2.5)
       end
       hold off 

       title();
       xlabel('Time(min)');
       ylabel('Glucocorticoid Inside The Cell(Molecules)');
   end

    if Plot_R
       figure,
       for i= 1:numCurves
           stairs(orgTime,newR(i,:))
           hold on
       end

       if PlotMean
           stairs(orgTime,RMean,'r','LineWidth',2.5)
       end
       hold off 

       title();
       xlabel('Time(min)');
       ylabel('Receptor(Molecules)');
   end           

    if Plot_CC
       figure,
       for i= 1:numCurves
           stairs(orgTime,newCC(i,:))
           hold on
       end

       if PlotMean
           stairs(orgTime,CCMean,'r','LineWidth',2.5)
       end
       hold off 

       title();
       xlabel('Time(min)');
       ylabel('Receptor-Glucocorticoid Complex(Molecules)');
   end           

   if Plot_V
       figure,
       for i= 1:numCurves
           stairs(orgTime,newV(i,:))
           hold on
       end

       if PlotMean
           stairs(orgTime,VMean,'r','LineWidth',2.5)
       end
       hold off 

       title();
       xlabel('Time(min)');
       ylabel('Signal Production(Molecules)');
   end







Nickel removal system

Deterministic model

Equations

function y = EqNick(x,t)

%--Parameters---%

global gammaN Kp beta Kd Kx alfaR deltaR Kt deltaT n alfaP deltaP


No=x(1); % Niquel outside the cell 

%------ Variables%------%

Ni= x(2); %Nickel Inside the cell R=x(3); %RcnR (Repressor) T=x(4); %RcnR Tetramer P=x(5); %Porine 


                                %---Equations---%

dNo=-gammaN*(No-Ni) - Kp*P*No; dNi=gammaN*(No-Ni) + Kp*P*No- beta/(1 + (T/(Kd*(1+(Ni/Kx))^n))); dR=alfaR-deltaR*R - Kt*R^4; dT= Kt*R^4 - deltaT*T - beta/(1 + (T/(Kd*(1+(Ni/Kx))^n))); dP=alfaP - deltaP*P + beta/(1 + (T/(Kd*(1+(Ni/Kx))^n)));

y1(1)=dNo; y1(2)=dNi; y1(3)=dR; y1(4)=dT; y1(5)=dP;

y= y1';

end

Equation solver

%
%
clear all
clc
%
%
%
%---------Parameters------%
%
%
%
global gammaN Kp beta Kd Kx alfaR deltaR Kt deltaT n alfaP deltaP
%
gammaN=0.5034e-4; %Diffussion rate of Nickel (1/min)
%
Kp=0.000634; %Dynamic constant for the entrance of nickel to the cell
%
beta=0.166; %Porine maximum expression rate (nM/min)
%
Kd=276e-3; %Association constant DNA and repressor (nM)
%
Kx=25e-3; %Association constant of the repressor with nickel (nM)
%
alfaR= 5; %Represor basal production rate (nM/min)
%
deltaR=1/1200; % Represor destruction rate (1/min)
%
Kt=820e-3;% Rate constant for the formation of the tetramer (Giraldo et al)
%
deltaT=1/1200; %Tetramer destruction rate (!/min)
%
n=1; %Hill coefficient (cooperation constant)
%
alfaP=0.031; %Porine basal production rate (nM/min)
%
deltaP=1/1200;%Porine destruction rate (1/min)
%
%
%
yo=[0 0 0 0];
%
con=fsolve(@(x)CondIniciales(x),yo, optimset('display','iter','MaxIter',1000000,'algorithm','levenberg-marquardt','tolfun',1e-9));
%
assignin('base','ini',con);
%
%
cond=abs(con);
%
h=30; %Tiempo maximo
%
m=0.01; %Paso
l=(0:m:h);%Vector de tiemp
%
condI=[9.88e3 cond];
x=zeros(length(l),length(condI));
x(1,:)=condI;
%
for k=1:length(l)-1
%    
   xk=x(k,:); %Captura de la ultima posicion de la matirz, es decir, los valores actuales de las variables
%   
   k1=EqNick(xk,l(k)); %Primera pendiente del metodo de RK4
   k2=EqNick(xk+(m/2*k1)',l(k)+m/2); %Segunda pendiente del metodo de RK4
   k3=EqNick(xk+(m/2*k2)',l(k)+m/2); %Tercera pendiente del metodo de RK4
   k4=EqNick(xk+(m*k3)',l(k)+m); %Cuarta pendiente del metodo de RK4
%  
   xk1=xk+m/6*(k1+2*k2+2*k3+k4)'; %Calculo de nuevos valores para las
   %variables
%  
   %xk1=xk+m*ecuaDif(l(k),xk)'; %Method of Newton
%   
   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
%
No=x(:,1);
Ni=x(:,2);
assignin('base','No',No);
assignin('base','Nistable',Ni(length(Ni)));
disp(Ni(length(Ni)))
R=x(:,3);
assignin('base','Rstable',R(length(R)));
T=x(:,4);
assignin('base','Tstable',T(length(T)));
P=x(:,5);
assignin('base','Pstable',P(length(P)));
cond=[R(length(R)) T(length(T)) P(length(P))];
assignin('base','cond',cond);
figure(1)
plot(l,No,l,P)
legend('No','P')
xlabel('Time (min)')
ylabel('Concentration (nM)')
figure(2)
plot(l,Ni)
legend('Ni')
xlabel('Time (min)')
ylabel('Concentration (nM)')
figure(3)
plot(l,R)
legend('R')
xlabel('Time (min)')
ylabel('Concentration (nM)')
figure(4)
plot(l,T)
legend('T')
xlabel('Time (min)')
ylabel('Concentration (nM)')
figure(5)
plot(l,P)
legend('P')
xlabel('Time (min)')
ylabel('Concentration (nM)') 
%

Stochastic

Retrieved from "http://2013.igem.org/Team:Colombia_Uniandes/Scripting"