Team:Colombia Uniandes/Scripting

From 2013.igem.org

(Difference between revisions)
(→Stochastic)
 
(18 intermediate revisions not shown)
Line 1: Line 1:
{{Http://2013.igem.org/Team:Colombia Uniandes/Bootstrap}}
{{Http://2013.igem.org/Team:Colombia Uniandes/Bootstrap}}
-
 
-
 
<html>
<html>
-
<div class="container-fluid">
+
  <div class="container">
-
  <div class="row-fluid">
+
    <div class="span11">
-
    <div class="span3">
+
<center><h2>Scripting</h2></center>
-
<div class="alert alert-info">
+
<h4 id="chimi"> Glucocorticoid Detection System Deterministic model</h4>
-
<ul class="nav nav-list">
+
<div class="accordion" id="accordion2">
-
  <li class="nav-header">Index</li>
+
  <div class="accordion-group">
-
  <li class="active"><a href="https://2013.igem.org/Team:Colombia_Uniandes/Scripting">Scripting</a></li>
+
    <div class="accordion-heading">
-
<ul>
+
      <a class="accordion-toggle" data-toggle="collapse" data-parent="#accordion2" href="#collapseOne">
-
<li><a href="https://2013.igem.org/Team:Colombia_Uniandes/Scripting#chimi">Glucocorticoid detection system</a></li>
+
Equations
-
  <li><a href="https://2013.igem.org/Team:Colombia_Uniandes/Scripting#nico">Nickel removal system</a></li>
+
      </a>
-
</ul>
+
    </div>
-
</ul>
+
    <div id="collapseOne" class="accordion-body collapse">
 +
      <div class="accordion-inner">
-
</div>
+
<h4>Equations</h4>
-
  </div>
+
-
    <div class="span9">
+
-
<h2>Scripting</h2>
+
</html>
</html>
-
 
-
<h3 id="chimi"> Glucocorticoid Detection System </h3>
 
-
 
-
 
-
===Deterministic model===
 
-
 
-
====Equations====
 
-
 
  function y = EcuacionesGluco(t,x)
  function y = EcuacionesGluco(t,x)
  global gammaGR mGRIR mCC deltaGRI alfaR deltaR deltaCC betaCC k n deltaS H
  global gammaGR mGRIR mCC deltaGRI alfaR deltaR deltaCC betaCC k n deltaS H
Line 60: Line 48:
  %
  %
  end
  end
-
 
+
<html>
-
====Equation solver====
+
      </div>
-
 
+
    </div>
 +
  </div>
 +
  <div class="accordion-group">
 +
    <div class="accordion-heading">
 +
      <a class="accordion-toggle" data-toggle="collapse" data-parent="#accordion2" href="#collapseTwo">
 +
      Equation solver
 +
      </a>
 +
    </div>
 +
    <div id="collapseTwo" class="accordion-body collapse">
 +
      <div class="accordion-inner">
 +
<h4>Equation solver</h4>
 +
</html>
Line 68: Line 67:
  %
  %
  %
  %
-
  gammaGR= 0.1;   %Diffussion rate of glucocorticoid inside the cell (mm/min)
+
  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)
-
  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)
  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)
  deltaGRI=0.00833;  %Glucocorticoids Destruction rate inside the cell (1/min)
-
%
+
  alfaR= 0.8e3;     %Basal production rate of the receptor (umol/min)  
-
  alfaR= 0.8e3;           %Basal production rate of the receptor (umol/min)  
+
  deltaR=0.004166;   %Receptor destruction rate inside the cell (1/min)
-
%
+
-
  deltaR=0.004166;         %Receptor destruction rate inside the cell (1/min)
+
  deltaCC=0.004166;  % GRI-R complex Destruction rate (1/min)
  deltaCC=0.004166;  % GRI-R complex Destruction rate (1/min)
-
  betaCC=0.5e3;     % GRI-R complex maximum expression rate (umol/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)
  k=0.05e3;          %Hill's constant for the GRI-R complex dimmer binding to his respective region (umol)
  n=2;              %Hill coefficient (cooperation constant)
  n=2;              %Hill coefficient (cooperation constant)
-
  deltaS=0.04166;   %Signal destruction rate (1/min)
+
  deltaS=0.04166;   %Signal destruction rate (1/min)
  H=2;              %Correction constant for the signal
  H=2;              %Correction constant for the signal
  %
  %
  %
  %
  %
  %
-
  h=60; %Tiempo maximo
+
  h=60; %maximum time
-
%
+
  m=0.01; %step length [s]
-
  m=0.01; %Longitud de paso [s]
+
  t=0:m:h; %time vector
-
%
+
-
  t=0:m:h; %Vector tiempo
+
  %
  %
  xi=[0 0 0 0];
  xi=[0 0 0 0];
Line 111: Line 103:
  for u=1:length(l)-1
  for u=1:length(l)-1
  %     
  %     
-
     xk=x(u,:); %Captura de la ultima posicion de la matirz, es decir, los
+
     xk=x(u,:); %Takes the last position of the matrix (i.e. the actual values of the variables)
-
    %valores actuales de las variables
+
  %     
  %     
-
    k1=EcuacionesGluco(l(u),xk); %Primera pendiente del metodo de RK4
+
%  k1=EcuacionesGluco(l(u),xk); %first slope of the RK4 method
-
    k2=EcuacionesGluco(l(u)+m/2,xk+(m/2*k1)'); %Segunda pendiente del metodo de RK4
+
%  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)'); %Tercera pendiente del metodo de RK4
+
%  k3=EcuacionesGluco(l(u)+m/2,xk+(m/2*k2)'); %third slope of the RK4 method
-
    k4=EcuacionesGluco(l(u)+m,xk+(m*k3)'); %Cuarta pendiente del metodo de RK4
+
%  k4=EcuacionesGluco(l(u)+m,xk+(m*k3)'); %ffourth slope of the RK4 method
  %   
  %   
-
    xk1=xk+m/6*(k1+2*k2+2*k3+k4)'; %Calculo de nuevos valores para las
+
%  xk1=xk+m/6*(k1+2*k2+2*k3+k4)'; %new vectors for the variables
-
    %variables
+
  %   
  %   
  %       
  %       
-
    xk2=zeros(1,length(xk1));
+
%  xk2=zeros(1,length(xk1));
  %   
  %   
  %   
  %   
-
    for p=1:length(xk1)
+
%  for p=1:length(xk1)
  %       
  %       
-
        if(xk1(p)<0.00000001)
+
%      if(xk1(p)<0.00000001)
  %           
  %           
-
            xk2(p)=0;
+
%          xk2(p)=0;
-
        else
+
%      else
  %           
  %           
-
            xk2(p)=xk1(p);
+
%          xk2(p)=xk1(p);
-
        end
+
%      end
  %       
  %       
-
    end
+
%  end
  %   
  %   
  %   
  %   
-
    x(u+1,:)=xk2; %Actualizacion del nuevo vector de variables en la matriz
+
%  x(u+1,:)=xk2; %Actualization of the new vector of variables int the martrix
  %   
  %   
  %   
  %   
Line 149: Line 139:
  for j=1:length(l)
  for j=1:length(l)
  %   
  %   
-
    if (l(j)<(10) || l(j)>(30))
+
%  if (l(j)<(10) || l(j)>(30))
  %       
  %       
-
        GRO(j)=155;
+
%      GRO(j)=155;
  %       
  %       
-
    else
+
%  else
  %       
  %       
-
        GRO(j)=155*1.3;
+
%      GRO(j)=155*1.3;
  %       
  %       
  %       
  %       
-
    end
+
%  end
  %   
  %   
  %   
  %   
Line 187: Line 177:
  plot(l,GRI)%,l,GRO)
  plot(l,GRI)%,l,GRO)
  legend('GRI')%,'Glucocorticoid')
  legend('GRI')%,'Glucocorticoid')
 +
%
 +
<html>   
 +
</div>
 +
    </div>
 +
  </div>
 +
<div class="accordion-group">
 +
    <div class="accordion-heading">
 +
      <a class="accordion-toggle" data-toggle="collapse" data-parent="#accordion3" href="#collapseThree">
 +
    Stochastic
 +
      </a>
 +
    </div>
 +
    <div id="collapseThree" class="accordion-body collapse">
 +
      <div class="accordion-inner">
 +
<h4>Stochastic</h4>
 +
</html>
  %
  %
-
 
+
%Stochastic Simulation of the Chemical Reactions of the Glucocorticoid
-
===Stochastic===
+
%System
-
 
+
%
-
%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
-
 
+
%
-
    %%----Parameters of the System----
+
%%-----Initial Conditions of the glucocorticoid system-----
-
    global gammaGR mGRIR mCC deltaGRI alfaR deltaR deltaCC betaCC k n deltaS H
+
%Run CondIndGluco.m with the parameters of the system.
-
 
+
%
-
    %%-----Initial Conditions of the glucocorticoid system-----
+
con=condiciones();
-
    %Run CondIndGluco.m with the parameters of the system.
+
%
-
 
+
GRO_0= round(10000); %Glucocorticoid Outside the Cell
-
    con=condiciones();
+
GRI_0 = round(con(1)); %Glucocorticoid Inside the Cell
-
 
+
R_0= round(con(2)); %Receptor
-
    GRO_0= round(10000); %Glucocorticoid Outside the Cell
+
CC_0= round(con(3)); %Receptor- Glucocorticoid Complex
-
    GRI_0 = round(con(1)); %Glucocorticoid Inside the Cell
+
V_0= round(con(4));%Signal
-
    R_0= round(con(2)); %Receptor
+
%
-
    CC_0= round(con(3)); %Receptor- Glucocorticoid Complex
+
%%----Simulation Constants----
-
    V_0= round(con(4));%Signal
+
numCells = 10;
-
 
+
numEvents= 60000;
-
    %%----Simulation Constants----
+
numCurves = 5; % Restriction(numCurves<numCells)
-
    numCells = 10;
+
Step=0.1;
-
    numEvents= 60000;
+
tMax=45; %Minutes
-
    numCurves = 5; % Restriction(numCurves<numCells)
+
%
-
    Step=0.1;
+
%Plot the behaviour of the molecule 1/0 for yes/no
-
    tMax=45; %Minutes
+
Plot_GRO = 1;             
-
 
+
Plot_GRI=1;             
-
    %Plot the behaviour of the molecule 1/0 for yes/no
+
Plot_R =1;             
-
    Plot_GRO = 1;             
+
Plot_CC = 1;             
-
    Plot_GRI=1;             
+
Plot_V = 1;   
-
    Plot_R =1;             
+
PlotMean=1;           
-
    Plot_CC = 1;             
+
%
-
    Plot_V = 1;   
+
%%----Stochastic Simulation----
-
    PlotMean=1;           
+
%Variables of the System
-
 
+
GRO= zeros(numCells,numEvents);
-
    %%----Stochastic Simulation----
+
GRI= zeros(numCells,numEvents);
-
    %Variables of the System
+
R= zeros(numCells,numEvents);
-
    GRO= zeros(numCells,numEvents);
+
CC= zeros(numCells,numEvents);
-
    GRI= zeros(numCells,numEvents);
+
V = zeros(numCells,numEvents);
-
    R= zeros(numCells,numEvents);
+
%
-
    CC= zeros(numCells,numEvents);
+
t = zeros(numCells,numEvents); %Time
-
    V = zeros(numCells,numEvents);
+
tr = zeros(numCells,numEvents);
-
 
+
%Initial Conditions
-
    t = zeros(numCells,numEvents); %Time
+
GRO(:,1)=GRO_0*ones(numCells,1);
-
    tr = zeros(numCells,numEvents);
+
GRI(:,1)=GRI_0*ones(numCells,1);
-
    %Initial Conditions
+
R(:,1)=R_0*ones(numCells,1);
-
    GRO(:,1)=GRO_0*ones(numCells,1);
+
CC(:,1)=CC_0*ones(numCells,1);
-
    GRI(:,1)=GRI_0*ones(numCells,1);
+
V(:,1)=V_0*ones(numCells,1);
-
    R(:,1)=R_0*ones(numCells,1);
+
%
-
    CC(:,1)=CC_0*ones(numCells,1);
+
%Stochastic Simulation
-
    V(:,1)=V_0*ones(numCells,1);
+
%
-
 
+
for S_ns =1:numCells
-
    %Stochastic Simulation
+
%
-
 
+
%
-
    for S_ns =1:numCells
+
    for S_ne = 1:numEvents
-
 
+
        %Random numbers drawn from uniform distribution between zero and
-
 
+
        %one
-
        for S_ne = 1:numEvents
+
        rand1 = rand(1); %This determine time
-
            %Random numbers drawn from uniform distribution between zero and
+
        rand2 = rand(1); %This determine Event
-
            %one
+
%
-
            rand1 = rand(1); %This determine time
+
        %Events (Description in attached document)
-
            rand2 = rand(1); %This determine Event
+
        e1= gammaGR*GRO(S_ns,S_ne);
-
 
+
        e2= gammaGR*GRI(S_ns,S_ne);
-
            %Events (Description in attached document)
+
        e3= mGRIR*GRI(S_ns,S_ne)*R(S_ns,S_ne);
-
            e1= gammaGR*GRO(S_ns,S_ne);
+
        e4= mCC*CC(S_ns,S_ne);
-
            e2= gammaGR*GRI(S_ns,S_ne);
+
        e5= deltaGRI*GRI(S_ns,S_ne);
-
            e3= mGRIR*GRI(S_ns,S_ne)*R(S_ns,S_ne);
+
        e6= alfaR;
-
            e4= mCC*CC(S_ns,S_ne);
+
        e7= deltaR*R(S_ns,S_ne);
-
            e5= deltaGRI*GRI(S_ns,S_ne);
+
        e8= deltaCC*CC(S_ns,S_ne);
-
            e6= alfaR;
+
        %e8= (betaCC*((CC(S_ns,S_ne))^n))/(k^n)+ ((CC(S_ns,S_ne))^n);
-
            e7= deltaR*R(S_ns,S_ne);
+
        e9= H*((betaCC*((CC(S_ns,S_ne))^n))/((k^n)+ ((CC(S_ns,S_ne))^n)));
-
            e8= deltaCC*CC(S_ns,S_ne);
+
        e10= deltaS*V(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)));
+
        %Vector of Events
-
            e10= deltaS*V(S_ns,S_ne);
+
        vecEvents = [e1,e2,e3,e4,e5,e6,e7,e8,e9,e10];
-
 
+
%
-
            %Vector of Events
+
        %Sum of Events
-
            vecEvents = [e1,e2,e3,e4,e5,e6,e7,e8,e9,e10];
+
        sumEvents = sum(vecEvents);
-
 
+
%
-
            %Sum of Events
+
        %Waitiing time until the next ocurrence of a reaction time
-
            sumEvents = sum(vecEvents);
+
        %Probabilistic Distribution
-
 
+
        time= (1/sumEvents)*log(1/rand1);
-
            %Waitiing time until the next ocurrence of a reaction time
+
        tr(S_ns,S_ne)=time;
-
            %Probabilistic Distribution
+
        t(S_ns,S_ne+1)= t(S_ns,S_ne)+ time;
-
            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));
-
            %Glucocorticoid outside the cell         
+
        if(t(S_ns,S_ne)<130)
-
 
+
            GRI(S_ns,S_ne)=funcImpulso(t(S_ns,S_ne+1))*957/2575;
-
            GRO(S_ns,S_ne+1)=funcImpulso(t(S_ns,S_ne+1));
+
        end
-
            if(t(S_ns,S_ne)<130)
+
%
-
                GRI(S_ns,S_ne)=funcImpulso(t(S_ns,S_ne+1))*957/2575;
+
        %Simulation Creation/Destruction of each molecule
-
            end
+
        if rand2<(vecEvents(1)/sumEvents)
-
 
+
%
-
            %Simulation Creation/Destruction of each molecule
+
            GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne)+1;
-
            if rand2<(vecEvents(1)/sumEvents)
+
            R(S_ns, S_ne+1)= R(S_ns,S_ne);
-
 
+
            CC(S_ns, S_ne+1)= CC(S_ns,S_ne);
-
                GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne)+1;
+
            V(S_ns, S_ne+1) = V(S_ns,S_ne);
-
                R(S_ns, S_ne+1)= R(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);
                 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);
                 R(S_ns, S_ne+1)= R(S_ns,S_ne);
-
                 CC(S_ns, S_ne+1)= CC(S_ns,S_ne);
+
                GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne);
-
                 V(S_ns, S_ne+1) = V(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
 +
%
 +
%
 +
%
 +
%
 +
%
 +
%
 +
%
 +
%
 +
%
 +
%
 +
%
 +
%
 +
%
 +
<html>
 +
      </div>
 +
    </div>
 +
  </div>
 +
</div>
-
            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
+
<h4 id="nico"> Nickel removal system Deterministic model</h4>
-
                if CC(S_ns,S_ne)<=50
+
<div class="accordion" id="accordion3">
-
                  CC(S_ns,S_ne+1)=0;
+
  <div class="accordion-group">
-
                  GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne);
+
    <div class="accordion-heading">
-
                  R(S_ns, S_ne+1)= R(S_ns,S_ne);  
+
      <a class="accordion-toggle" data-toggle="collapse" data-parent="#accordion1" href="#collapse1">
-
                else
+
Equations
-
                  CC(S_ns, S_ne+1)= CC(S_ns,S_ne)-50;
+
      </a>
-
                  GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne)+50;
+
    </div>
-
                  R(S_ns, S_ne+1)= R(S_ns,S_ne)+50;  
+
    <div id="collapse1" class="accordion-body collapse">
-
                end
+
      <div class="accordion-inner">
 +
</html>
 +
<h4>Equations</h4>
 +
%
 +
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
-
                V(S_ns, S_ne+1) = V(S_ns,S_ne);
+
<html>
 +
      </div>
 +
    </div>
 +
  </div>
 +
  <div class="accordion-group">
 +
    <div class="accordion-heading">
 +
      <a class="accordion-toggle" data-toggle="collapse" data-parent="#accordion2" href="#collapse2">
 +
        Equation solver
 +
      </a>
 +
    </div>
 +
    <div id="collapse2" class="accordion-body collapse">
 +
      <div class="accordion-inner">
 +
</html>
 +
<h3>Equation solver</h3>
 +
%
 +
%
 +
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)')
 +
%
-
            elseif rand2< (sum(vecEvents(1:5)))/sumEvents
+
<html>
 +
      </div>
 +
    </div>
 +
  </div>
 +
<div class="accordion-group">
 +
    <div class="accordion-heading">
 +
      <a class="accordion-toggle" data-toggle="collapse" data-parent="#accordion3" href="#collapse3">
 +
Stochastic
 +
      </a>
 +
    </div>
 +
    <div id="collapse3" class="accordion-body collapse">
 +
      <div class="accordion-inner">
 +
<h3>Stochastic</h3>
-
                if GRI(S_ns,S_ne)<=0
+
</html>
-
                  GRI(S_ns,S_ne+1)=0;
+
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +
    % 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
                 else
-
                  GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne)-10;
+
                    T(s_NS,s_NE) = 0;
                 end
                 end
-
 
+
                 P(s_NS,s_NE+1)= P(s_NS,s_NE)+changeCreation;
-
 
+
      %
-
                 R(S_ns, S_ne+1)= R(S_ns,S_ne);
+
             elseif s2 < (sum(VecEvents(1:2)/SumEvents))
-
                CC(S_ns, S_ne+1)= CC(S_ns,S_ne);
+
      %
-
                V(S_ns, S_ne+1) = V(S_ns,S_ne);  
+
                 R(s_NS,s_NE+1)= R(s_NS,s_NE)+changeCreation;
-
 
+
                 T(s_NS,s_NE+1)= T(s_NS,s_NE);
-
             elseif rand2< (sum(vecEvents(1:6)))/sumEvents
+
                 P(s_NS,s_NE+1)= P(s_NS,s_NE);
-
 
+
    %
-
                GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne);
+
             elseif s2 < (sum(VecEvents(1:3)/SumEvents))
-
                 R(S_ns, S_ne+1)= R(S_ns,S_ne)+600;
+
      %
-
                 CC(S_ns, S_ne+1)= CC(S_ns,S_ne);
+
                if R(s_NS,s_NE)>0
-
                 V(S_ns, S_ne+1) = V(S_ns,S_ne);
+
                    R(s_NS,s_NE+1)= R(s_NS,s_NE)-changeDecrease;
-
 
+
-
             elseif rand2< (sum(vecEvents(1:7)))/sumEvents
+
-
 
+
-
                if R(S_ns,S_ne)<=0
+
-
                      R(S_ns,S_ne+1)=0;
+
                 else
                 else
-
                     R(S_ns, S_ne+1)= R(S_ns,S_ne)-300;
+
                     R(s_NS,s_NE) = 0;
                 end
                 end
-
 
+
                 T(s_NS,s_NE+1)= T(s_NS,s_NE);
-
 
+
                 P(s_NS,s_NE+1)= P(s_NS,s_NE);
-
                 GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne);
+
    %
-
                 CC(S_ns, S_ne+1)= CC(S_ns,S_ne);
+
             elseif s2 < (sum(VecEvents(1:4)/SumEvents))
-
                V(S_ns, S_ne+1) = V(S_ns,S_ne);
+
      %
-
 
+
              if R(s_NS,s_NE)>0
-
             elseif rand2< (sum(vecEvents(1:8)))/sumEvents
+
                    R(s_NS,s_NE+1)= R(s_NS,s_NE)-changeDecrease;
-
                if CC(S_ns,S_ne)<=35
+
-
                  CC(S_ns,S_ne+1)=round(CC(S_ns,S_ne)*0.75);
+
                 else
                 else
-
                  CC(S_ns, S_ne+1)= CC(S_ns,S_ne)-35;
+
                    R(s_NS,s_NE) = 0;
-
 
+
                 end
                 end
-
                 GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne);
+
                 T(s_NS,s_NE+1)= T(s_NS,s_NE)+changeCreation;
-
                R(S_ns, S_ne+1)= R(S_ns,S_ne);                  
+
                 P(s_NS,s_NE+1)= P(s_NS,s_NE);
-
                 V(S_ns, S_ne+1) = V(S_ns,S_ne);
+
      %
-
 
+
             elseif s2 < (sum(VecEvents(1:5)/SumEvents))
-
             elseif rand2< (sum(vecEvents(1:9)))/sumEvents
+
      %
-
 
+
                 R(s_NS,s_NE+1)= R(s_NS,s_NE);
-
                GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne);
+
                 if T(s_NS,s_NE)>0
-
                 R(S_ns, S_ne+1)= R(S_ns,S_ne);
+
                    T(s_NS,s_NE+1)= T(s_NS,s_NE)-changeDecrease;
-
                 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
                 else
-
                  V(S_ns, S_ne+1)= V(S_ns,S_ne)-500;
+
                    T(s_NS,s_NE) = 0;
                 end
                 end
-
 
+
                P(s_NS,s_NE+1)= P(s_NS,s_NE);
-
                 GRI(S_ns,S_ne+1)= GRI(S_ns,S_ne);
+
      %
-
                 R(S_ns, S_ne+1)= R(S_ns,S_ne);
+
            elseif s2 < (sum(VecEvents(1:6)/SumEvents))
-
                 CC(S_ns, S_ne+1)= CC(S_ns,S_ne);
+
      %
-
 
+
                 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
         end
     end
     end
-
 
+
    %
 +
    % Nickel Inside the Cell
 +
    %
 +
    Ni = NickPerPorine*P;
 +
  %
 +
    %
     %% ------------- Time Regularization -----------------------------
     %% ------------- Time Regularization -----------------------------
-
     time=t;
+
     time = t;
-
 
+
     %
-
     [orgTime, newGRO] = regIntervalFixed2(time, GRO, Step, tMax);
+
     [~, newNi] = regIntervalFixed2(time, Ni, Step, TimeMax);
-
     [orgTime, newGRI] = regIntervalFixed2(time,GRI, Step, tMax);
+
     [~, newR] = regIntervalFixed2(time, R, Step, TimeMax);
-
     [orgTime, newR] = regIntervalFixed2(time, R, Step, tMax);
+
     [~, newT] = regIntervalFixed2(time, T, Step, TimeMax);
-
     [orgTime, newCC] = regIntervalFixed2(time, CC, Step, tMax);
+
     [newTime, newP] = regIntervalFixed2(time, P, Step, TimeMax);
-
     [orgTime, newV] = regIntervalFixed2(time, V, Step, tMax);
+
    %
-
 
+
     %
-
     GROMean=mean(newGRO);
+
     NiMean=mean(newNi);
-
     GRIMean=mean(newGRI);
+
     RMean=mean(newR);
     RMean=mean(newR);
-
     CCMean=mean(newCC);
+
     TMean=mean(newT);
-
     VMean=mean(newV);
+
     %
-
 
+
     %% ----------------------- Plotting ----------------------------------
-
 
+
     %
-
 
+
     %
-
     %% Plotting  
+
     if Plot_Ni
-
     % 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,
         figure,
-
 
+
         for s_NP = 1:NumbCurvesPlotted
-
         for i= 1:numCurves
+
             stairs(newTime,newNi(s_NP,:))
-
             stairs(orgTime,newGRO(i,:))
+
             hold on
             hold on
         end
         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
         if PlotMean
-
             stairs(orgTime,GRIMean,'r','LineWidth',2.5)
+
             stairs(newTime,NiMean,'r');
         end
         end
-
         hold off  
+
         hold off
-
 
+
         xlabel('Time (min)');
-
        title('');
+
         ylabel('Nickel Inside the cell (Molecules)');
-
         xlabel('Time(min)');
+
-
         ylabel('Glucocorticoid Inside The Cell(Molecules)');
+
     end
     end
-
 
+
  %
-
    if Plot_R
+
    if Plot_R
         figure,
         figure,
-
         for i= 1:numCurves
+
         for s_NP = 1:NumbCurvesPlotted
-
             stairs(orgTime,newR(i,:))
+
             stairs(newTime,newR(s_NP,:))
             hold on
             hold on
         end
         end
-
 
         if PlotMean
         if PlotMean
-
             stairs(orgTime,RMean,'r','LineWidth',2.5)
+
             stairs(newTime,mean(newR),'r');
         end
         end
-
         hold off  
+
         hold off
-
 
+
         xlabel('Time (min)');
-
        title('');
+
         ylabel('Repressor (Molecules)');
-
         xlabel('Time(min)');
+
     end
-
         ylabel('Receptor(Molecules)');
+
  %
-
     end          
+
    if Plot_T
-
 
+
-
    if Plot_CC
+
         figure,
         figure,
-
         for i= 1:numCurves
+
         for s_NP = 1:NumbCurvesPlotted
-
             stairs(orgTime,newCC(i,:))
+
             stairs(newTime,newT(s_NP,:))
             hold on
             hold on
         end
         end
-
 
         if PlotMean
         if PlotMean
-
             stairs(orgTime,CCMean,'r','LineWidth',2.5)
+
             stairs(newTime,mean(newT),'r');
         end
         end
-
         hold off  
+
         hold off
-
 
+
         xlabel('Time (min)');
-
        title('');
+
         ylabel('Tetramer (Molecules)');
-
         xlabel('Time(min)');
+
     end
-
         ylabel('Receptor-Glucocorticoid Complex(Molecules)');
+
    %
-
     end          
+
     if Plot_P
-
 
+
-
     if Plot_V
+
         figure,
         figure,
-
         for i= 1:numCurves
+
         for s_NP = 1:NumbCurvesPlotted
-
             stairs(orgTime,newV(i,:))
+
             stairs(newTime,newP(s_NP,:))
             hold on
             hold on
         end
         end
-
 
         if PlotMean
         if PlotMean
-
             stairs(orgTime,VMean,'r','LineWidth',2.5)
+
             stairs(newTime,mean(newP),'r');
         end
         end
-
         hold off  
+
         hold off
-
 
+
         xlabel('Time (min)');
-
        title('');
+
         ylabel('Porines (Molecules)');
-
         xlabel('Time(min)');
+
     end
-
         ylabel('Signal Production(Molecules)');
+
-
     end  
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
 +
<html>
-
<h3 id="nico"> Nickel removal system </h3>
 
-
===Deterministic model===
 
-
====Equations====
+
      </div>
-
 
+
    </div>
-
function y = EqNick(x,t)
+
  </div>
-
 
+
</div>
-
%--Parameters---%
+
  </div>
-
 
+
</div>
-
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===
+

Latest revision as of 21:38, 26 September 2013

Scripting

Glucocorticoid Detection System Deterministic model

Equations

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

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

%
%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

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

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

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % 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


Retrieved from "http://2013.igem.org/Team:Colombia_Uniandes/Scripting"