Team:Colombia Uniandes/Scripting

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