Team:Colombia Uniandes/Scripting

Scripting

Contents 1 Glucocorticoid Detection System 2 Deterministic model 2.1 Equations 2.2 Equation solver 3 Stochastic 4 Nickel removal system 5 Deterministic model 5.1 Equations 5.2 Equation solver 6 Stochastic

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