"
Team:Colombia Uniandes/Scripting
From 2013.igem.org
Revision as of 22:21, 25 September 2013 by D.olivera1320 (Talk | contribs)
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)')
