"
Team:Colombia Uniandes/Scripting
From 2013.igem.org
Revision as of 21:38, 26 September 2013 by D.olivera1320 (Talk | contribs)
Scripting
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; %maximum time
m=0.01; %step length [s]
t=0:m:h; %time vector
%
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,:); %Takes the last position of the matrix (i.e. the actual values of the variables)
%
% k1=EcuacionesGluco(l(u),xk); %first slope of the RK4 method
% k2=EcuacionesGluco(l(u)+m/2,xk+(m/2*k1)'); %second slope of the RK4 method
% k3=EcuacionesGluco(l(u)+m/2,xk+(m/2*k2)'); %third slope of the RK4 method
% k4=EcuacionesGluco(l(u)+m,xk+(m*k3)'); %ffourth slope of the RK4 method
%
% xk1=xk+m/6*(k1+2*k2+2*k3+k4)'; %new vectors for the 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; %Actualization of the new vector of variables int the martrix
%
%
%
%
%
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
%
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 Simulation of the Chemical Reaction of the Nickel Absortion
% System
%
%
% Author: Igem Team Colombia 2013.
%
%% ----------------------- Parameters --------------------------------
%
clear all;close all;clc
%
global gammaN Kp beta Kd Kx alfaR deltaR Kt deltaT n alfaP deltaP
%
gammaN=0.5034; %Diffussion rate of Nickel (1/min)
%
Kp=0.00004634; %Dynamic constant for the entrance of nickel to the cell
%
beta=1.00e-1; %Porine maximum expression rate (nM/min)
%
Kd=1.66e2; %Association constant DNA and repressor (nM)
%
Kx=1.51e1; %Association constant of the repressor with nickel (nM)
%
alfaR= 6.02e-2; %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
%
deltaT=1/1200; %Tetramer destruction rate (!/min)
%
n=1; %Hill coefficient (cooperation constant)
%
alfaP=1.87e-2; %Porine basal production rate (nM/min)
%
deltaP=1/1200;%Porine destruction rate (1/min)
%
%% ---------------------- Initial Conditions ----------------------
%
% Taken from running CondIniciales.m with the same parameters above.
con=CondInicialesSto();
condI=[0 round(con)];
%
%condI=[100 0 0.590218512408079,118.952234340367,37.6611256746707];
%
Ni_0 = condI(1);
R_0 = condI(2);
T_0 = condI(3);
P_0 = condI(4);
%
%
%% --------------------- Simulation Constants ----------------------
%
NumbSimul = 100; % Number of Different Simulation
NumbEvent = 100000; % Number of Total Events happening along the simulation.
NumbCurvesPlotted = 7; % Number of Curves shown on each plot.(NumbCurvesPlotted < NumbSimul)
Step = 0.5; % Step between time intervals
TimeMax = 100; % Maximum time of the simulation.
PlotMean = 1; % Plot the curve of the mean of each simulation? 1/0 for yes/no
%
%
Plot_Ni = 1; % Plot the graph of Ni? 1/0 for yes/no
Plot_R = 1; % Plot the graph of R? 1/0 for yes/no
Plot_T = 1; % Plot the graph of T? 1/0 for yes/no
Plot_P = 1; % Plot the graph of P? 1/0 for yes/no
%
changeCreation = 1000; % Number of creation of proteins in certain events
changeDecrease = 1; % Number of molecules decreased in certain events
%
% Number of Nickel molecules inside dthe cell per number of Porine
% Molecules
NickPerPorine = 1/(7*(1/1e12)*(6.02e23)*(1.66e-21)*(40.1*1000));
%% --------------------- Stochastic Simulation ---------------------
%
% Variables
%
Ni = zeros(NumbSimul,NumbEvent);
R = zeros(NumbSimul,NumbEvent);
T = zeros(NumbSimul,NumbEvent);
P = zeros(NumbSimul,NumbEvent);
% Initial Conditions
%
Ni(:,1)=Ni_0*ones(NumbSimul,1);
R(:,1)=R_0*ones(NumbSimul,1);
T(:,1)=T_0*ones(NumbSimul,1);
P(:,1)=P_0*ones(NumbSimul,1);
%
% Time
t = zeros(NumbSimul,NumbEvent);
%
% Simulation
%
for s_NS = 1:NumbSimul
for s_NE = 1:NumbEvent
s1 = rand; % Random Number between 0 and 1 for time choosing
s2 = rand; % Random Number between 0 and 1 for event choosing
%
% Events (See document attached for better explanation)
%
Ev1 = beta/(1+(T(s_NS,s_NE)/(Kd*(1+(Ni(s_NS,s_NE)/Kx)^n))));
Ev2 = alfaR;
Ev3 = deltaR*R(s_NS,s_NE);
Ev4 = Kt*R(s_NS,s_NE)^4;
Ev5 = deltaT*T(s_NS,s_NE);
Ev6 = alfaP;
Ev7 = deltaP*P(s_NS,s_NE);
%
% Vector of Events
VecEvents = [Ev1,Ev2,Ev3,Ev4,Ev5,Ev6,Ev7];
%
SumEvents = sum(VecEvents);
%
% Chosen time at random based on the probability distribution
Time = (1/(SumEvents))*log(1/s1);
%
t(s_NS,s_NE+1) = t(s_NS,s_NE)+ Time;
%
% Event Choosing at random based on each event's weight. What
% happens on each event is explained better in the attached doc.
%
if s2 < (sum(VecEvents(1:1)/SumEvents))
%
R(s_NS,s_NE+1)= R(s_NS,s_NE);
if T(s_NS,s_NE)>0
T(s_NS,s_NE+1)= T(s_NS,s_NE)-changeDecrease;
else
T(s_NS,s_NE) = 0;
end
P(s_NS,s_NE+1)= P(s_NS,s_NE)+changeCreation;
%
elseif s2 < (sum(VecEvents(1:2)/SumEvents))
%
R(s_NS,s_NE+1)= R(s_NS,s_NE)+changeCreation;
T(s_NS,s_NE+1)= T(s_NS,s_NE);
P(s_NS,s_NE+1)= P(s_NS,s_NE);
%
elseif s2 < (sum(VecEvents(1:3)/SumEvents))
%
if R(s_NS,s_NE)>0
R(s_NS,s_NE+1)= R(s_NS,s_NE)-changeDecrease;
else
R(s_NS,s_NE) = 0;
end
T(s_NS,s_NE+1)= T(s_NS,s_NE);
P(s_NS,s_NE+1)= P(s_NS,s_NE);
%
elseif s2 < (sum(VecEvents(1:4)/SumEvents))
%
if R(s_NS,s_NE)>0
R(s_NS,s_NE+1)= R(s_NS,s_NE)-changeDecrease;
else
R(s_NS,s_NE) = 0;
end
T(s_NS,s_NE+1)= T(s_NS,s_NE)+changeCreation;
P(s_NS,s_NE+1)= P(s_NS,s_NE);
%
elseif s2 < (sum(VecEvents(1:5)/SumEvents))
%
R(s_NS,s_NE+1)= R(s_NS,s_NE);
if T(s_NS,s_NE)>0
T(s_NS,s_NE+1)= T(s_NS,s_NE)-changeDecrease;
else
T(s_NS,s_NE) = 0;
end
P(s_NS,s_NE+1)= P(s_NS,s_NE);
%
elseif s2 < (sum(VecEvents(1:6)/SumEvents))
%
R(s_NS,s_NE+1)= R(s_NS,s_NE);
T(s_NS,s_NE+1)= T(s_NS,s_NE);
P(s_NS,s_NE+1)= P(s_NS,s_NE)+changeCreation;
%
elseif s2 < (sum(VecEvents(1:7)/SumEvents))
%
R(s_NS,s_NE+1)= R(s_NS,s_NE);
T(s_NS,s_NE+1)= T(s_NS,s_NE);
if P(s_NS,s_NE)>0
P(s_NS,s_NE+1)= P(s_NS,s_NE)-changeDecrease;
else
P(s_NS,s_NE) = 0;
end;
end
%
end
end
%
% Nickel Inside the Cell
%
Ni = NickPerPorine*P;
%
%
%% ------------- Time Regularization -----------------------------
time = t;
%
[~, newNi] = regIntervalFixed2(time, Ni, Step, TimeMax);
[~, newR] = regIntervalFixed2(time, R, Step, TimeMax);
[~, newT] = regIntervalFixed2(time, T, Step, TimeMax);
[newTime, newP] = regIntervalFixed2(time, P, Step, TimeMax);
%
%
NiMean=mean(newNi);
RMean=mean(newR);
TMean=mean(newT);
%
%% ----------------------- Plotting ----------------------------------
%
%
if Plot_Ni
figure,
for s_NP = 1:NumbCurvesPlotted
stairs(newTime,newNi(s_NP,:))
hold on
end
if PlotMean
stairs(newTime,NiMean,'r');
end
hold off
xlabel('Time (min)');
ylabel('Nickel Inside the cell (Molecules)');
end
%
if Plot_R
figure,
for s_NP = 1:NumbCurvesPlotted
stairs(newTime,newR(s_NP,:))
hold on
end
if PlotMean
stairs(newTime,mean(newR),'r');
end
hold off
xlabel('Time (min)');
ylabel('Repressor (Molecules)');
end
%
if Plot_T
figure,
for s_NP = 1:NumbCurvesPlotted
stairs(newTime,newT(s_NP,:))
hold on
end
if PlotMean
stairs(newTime,mean(newT),'r');
end
hold off
xlabel('Time (min)');
ylabel('Tetramer (Molecules)');
end
%
if Plot_P
figure,
for s_NP = 1:NumbCurvesPlotted
stairs(newTime,newP(s_NP,:))
hold on
end
if PlotMean
stairs(newTime,mean(newP),'r');
end
hold off
xlabel('Time (min)');
ylabel('Porines (Molecules)');
end
Scripting
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; %maximum time
m=0.01; %step length [s]
t=0:m:h; %time vector
%
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,:); %Takes the last position of the matrix (i.e. the actual values of the variables)
%
% k1=EcuacionesGluco(l(u),xk); %first slope of the RK4 method
% k2=EcuacionesGluco(l(u)+m/2,xk+(m/2*k1)'); %second slope of the RK4 method
% k3=EcuacionesGluco(l(u)+m/2,xk+(m/2*k2)'); %third slope of the RK4 method
% k4=EcuacionesGluco(l(u)+m,xk+(m*k3)'); %ffourth slope of the RK4 method
%
% xk1=xk+m/6*(k1+2*k2+2*k3+k4)'; %new vectors for the 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; %Actualization of the new vector of variables int the martrix
%
%
%
%
%
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
%
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 Simulation of the Chemical Reaction of the Nickel Absortion
% System
%
%
% Author: Igem Team Colombia 2013.
%
%% ----------------------- Parameters --------------------------------
%
clear all;close all;clc
%
global gammaN Kp beta Kd Kx alfaR deltaR Kt deltaT n alfaP deltaP
%
gammaN=0.5034; %Diffussion rate of Nickel (1/min)
%
Kp=0.00004634; %Dynamic constant for the entrance of nickel to the cell
%
beta=1.00e-1; %Porine maximum expression rate (nM/min)
%
Kd=1.66e2; %Association constant DNA and repressor (nM)
%
Kx=1.51e1; %Association constant of the repressor with nickel (nM)
%
alfaR= 6.02e-2; %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
%
deltaT=1/1200; %Tetramer destruction rate (!/min)
%
n=1; %Hill coefficient (cooperation constant)
%
alfaP=1.87e-2; %Porine basal production rate (nM/min)
%
deltaP=1/1200;%Porine destruction rate (1/min)
%
%% ---------------------- Initial Conditions ----------------------
%
% Taken from running CondIniciales.m with the same parameters above.
con=CondInicialesSto();
condI=[0 round(con)];
%
%condI=[100 0 0.590218512408079,118.952234340367,37.6611256746707];
%
Ni_0 = condI(1);
R_0 = condI(2);
T_0 = condI(3);
P_0 = condI(4);
%
%
%% --------------------- Simulation Constants ----------------------
%
NumbSimul = 100; % Number of Different Simulation
NumbEvent = 100000; % Number of Total Events happening along the simulation.
NumbCurvesPlotted = 7; % Number of Curves shown on each plot.(NumbCurvesPlotted < NumbSimul)
Step = 0.5; % Step between time intervals
TimeMax = 100; % Maximum time of the simulation.
PlotMean = 1; % Plot the curve of the mean of each simulation? 1/0 for yes/no
%
%
Plot_Ni = 1; % Plot the graph of Ni? 1/0 for yes/no
Plot_R = 1; % Plot the graph of R? 1/0 for yes/no
Plot_T = 1; % Plot the graph of T? 1/0 for yes/no
Plot_P = 1; % Plot the graph of P? 1/0 for yes/no
%
changeCreation = 1000; % Number of creation of proteins in certain events
changeDecrease = 1; % Number of molecules decreased in certain events
%
% Number of Nickel molecules inside dthe cell per number of Porine
% Molecules
NickPerPorine = 1/(7*(1/1e12)*(6.02e23)*(1.66e-21)*(40.1*1000));
%% --------------------- Stochastic Simulation ---------------------
%
% Variables
%
Ni = zeros(NumbSimul,NumbEvent);
R = zeros(NumbSimul,NumbEvent);
T = zeros(NumbSimul,NumbEvent);
P = zeros(NumbSimul,NumbEvent);
% Initial Conditions
%
Ni(:,1)=Ni_0*ones(NumbSimul,1);
R(:,1)=R_0*ones(NumbSimul,1);
T(:,1)=T_0*ones(NumbSimul,1);
P(:,1)=P_0*ones(NumbSimul,1);
%
% Time
t = zeros(NumbSimul,NumbEvent);
%
% Simulation
%
for s_NS = 1:NumbSimul
for s_NE = 1:NumbEvent
s1 = rand; % Random Number between 0 and 1 for time choosing
s2 = rand; % Random Number between 0 and 1 for event choosing
%
% Events (See document attached for better explanation)
%
Ev1 = beta/(1+(T(s_NS,s_NE)/(Kd*(1+(Ni(s_NS,s_NE)/Kx)^n))));
Ev2 = alfaR;
Ev3 = deltaR*R(s_NS,s_NE);
Ev4 = Kt*R(s_NS,s_NE)^4;
Ev5 = deltaT*T(s_NS,s_NE);
Ev6 = alfaP;
Ev7 = deltaP*P(s_NS,s_NE);
%
% Vector of Events
VecEvents = [Ev1,Ev2,Ev3,Ev4,Ev5,Ev6,Ev7];
%
SumEvents = sum(VecEvents);
%
% Chosen time at random based on the probability distribution
Time = (1/(SumEvents))*log(1/s1);
%
t(s_NS,s_NE+1) = t(s_NS,s_NE)+ Time;
%
% Event Choosing at random based on each event's weight. What
% happens on each event is explained better in the attached doc.
%
if s2 < (sum(VecEvents(1:1)/SumEvents))
%
R(s_NS,s_NE+1)= R(s_NS,s_NE);
if T(s_NS,s_NE)>0
T(s_NS,s_NE+1)= T(s_NS,s_NE)-changeDecrease;
else
T(s_NS,s_NE) = 0;
end
P(s_NS,s_NE+1)= P(s_NS,s_NE)+changeCreation;
%
elseif s2 < (sum(VecEvents(1:2)/SumEvents))
%
R(s_NS,s_NE+1)= R(s_NS,s_NE)+changeCreation;
T(s_NS,s_NE+1)= T(s_NS,s_NE);
P(s_NS,s_NE+1)= P(s_NS,s_NE);
%
elseif s2 < (sum(VecEvents(1:3)/SumEvents))
%
if R(s_NS,s_NE)>0
R(s_NS,s_NE+1)= R(s_NS,s_NE)-changeDecrease;
else
R(s_NS,s_NE) = 0;
end
T(s_NS,s_NE+1)= T(s_NS,s_NE);
P(s_NS,s_NE+1)= P(s_NS,s_NE);
%
elseif s2 < (sum(VecEvents(1:4)/SumEvents))
%
if R(s_NS,s_NE)>0
R(s_NS,s_NE+1)= R(s_NS,s_NE)-changeDecrease;
else
R(s_NS,s_NE) = 0;
end
T(s_NS,s_NE+1)= T(s_NS,s_NE)+changeCreation;
P(s_NS,s_NE+1)= P(s_NS,s_NE);
%
elseif s2 < (sum(VecEvents(1:5)/SumEvents))
%
R(s_NS,s_NE+1)= R(s_NS,s_NE);
if T(s_NS,s_NE)>0
T(s_NS,s_NE+1)= T(s_NS,s_NE)-changeDecrease;
else
T(s_NS,s_NE) = 0;
end
P(s_NS,s_NE+1)= P(s_NS,s_NE);
%
elseif s2 < (sum(VecEvents(1:6)/SumEvents))
%
R(s_NS,s_NE+1)= R(s_NS,s_NE);
T(s_NS,s_NE+1)= T(s_NS,s_NE);
P(s_NS,s_NE+1)= P(s_NS,s_NE)+changeCreation;
%
elseif s2 < (sum(VecEvents(1:7)/SumEvents))
%
R(s_NS,s_NE+1)= R(s_NS,s_NE);
T(s_NS,s_NE+1)= T(s_NS,s_NE);
if P(s_NS,s_NE)>0
P(s_NS,s_NE+1)= P(s_NS,s_NE)-changeDecrease;
else
P(s_NS,s_NE) = 0;
end;
end
%
end
end
%
% Nickel Inside the Cell
%
Ni = NickPerPorine*P;
%
%
%% ------------- Time Regularization -----------------------------
time = t;
%
[~, newNi] = regIntervalFixed2(time, Ni, Step, TimeMax);
[~, newR] = regIntervalFixed2(time, R, Step, TimeMax);
[~, newT] = regIntervalFixed2(time, T, Step, TimeMax);
[newTime, newP] = regIntervalFixed2(time, P, Step, TimeMax);
%
%
NiMean=mean(newNi);
RMean=mean(newR);
TMean=mean(newT);
%
%% ----------------------- Plotting ----------------------------------
%
%
if Plot_Ni
figure,
for s_NP = 1:NumbCurvesPlotted
stairs(newTime,newNi(s_NP,:))
hold on
end
if PlotMean
stairs(newTime,NiMean,'r');
end
hold off
xlabel('Time (min)');
ylabel('Nickel Inside the cell (Molecules)');
end
%
if Plot_R
figure,
for s_NP = 1:NumbCurvesPlotted
stairs(newTime,newR(s_NP,:))
hold on
end
if PlotMean
stairs(newTime,mean(newR),'r');
end
hold off
xlabel('Time (min)');
ylabel('Repressor (Molecules)');
end
%
if Plot_T
figure,
for s_NP = 1:NumbCurvesPlotted
stairs(newTime,newT(s_NP,:))
hold on
end
if PlotMean
stairs(newTime,mean(newT),'r');
end
hold off
xlabel('Time (min)');
ylabel('Tetramer (Molecules)');
end
%
if Plot_P
figure,
for s_NP = 1:NumbCurvesPlotted
stairs(newTime,newP(s_NP,:))
hold on
end
if PlotMean
stairs(newTime,mean(newP),'r');
end
hold off
xlabel('Time (min)');
ylabel('Porines (Molecules)');
end
