Team:Colombia Uniandes/Scripting

From 2013.igem.org

(Difference between revisions)
(Stochastic)
(Deterministic model)
Line 201: Line 201:
<h3 id="nico"> Nickel removal system </h3>
<h3 id="nico"> Nickel removal system </h3>
-
==Deterministic model==
+
===Deterministic model===
-
===Equations===
+
====Equations====
-
===Equation solver===
+
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==
==Stochastic==

Revision as of 22:14, 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

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