Team:Colombia/Scripting

{{Http://2014.igem.org/Team:Colombia}}
<html>
<div class="span11" style="text-align: justify;">
<br>
<p>
<br>
<center><b> <font size="24">
<br>
Scripting
</font><br>
<h1  class="curs1"><font size="4">
Feel free to expand and scroll through the text boxes in order to further examine the code.
</font></h1></b></center>
<br>
<b><font color="#8A0808" size="3" >Deterministic Model</font></b>
<br>
This code creates the differential equations governing the concentration dinamics of each protein in our model, finds the steady state solutions and then solves them using the numerical aproximation method Runge-Kutta<br>
<textarea readonly rows = "20">
function y=det()
%
%--------------------------------------------%
%                 Runge Kutta 4 order aproximation           %
%--------------------------------------------%
%
h=1200;                                                             % Maximum time
m=0.005;                                                           % Step lenght [s]
xi=[1,1,1,1,1,1,1,1,1];                                        % 
y=fsolve(@CondIni,xi);                                        % 
conInd=[1,1,1,1,1,1,1,1,0];                                 % 
l=(0:m:h)';                                                           % Time vector
%
x=zeros(length(l),length(conInd));                        % Variable's Matrix, the variables in each column and time in each row
C=zeros(1,length(l));                                            % zero vector to be filled
x(1,:)=conInd;                                                      % 
%
for k=1:length(l)-1
    xk=x(k,:);                                                         % Capture of the last value in the matrix, the actual variable values
    k1=difeq(l(k),xk);                                             % First slope of the RK4 method
    k2=difeq(l(k)+m/2,xk+(m/2*k1)');                   % Second slope of the RK4 method
    k3=difeq(l(k)+m/2,xk+(m/2*k2)');                   % Third slope of the RK4 method
    k4=difeq(l(k)+m,xk+(m*k3)');                          % Fourth slope of the RK4 method
    xk1=xk+m/6*(k1+2*k2+2*k3+k4)';                % New value calculation for the variables
    %xk1=xk+m*difeq(l(k),xk)';                              % Another aproximation, Newton's method (comment line)                     
    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;                                                    % Variable vector actualization in the matrix
end
%
for j=1:length(l)
    if (l(j)>(300) & l(j)<(1000))
        C(j)=2.25;
    else
        C(j)=0;
    end
end
%

CS=x(:,1);                                                              % 
CSa=x(:,2);                                                            % 
Uf=x(:,3);                                                              % 
U=x(:,4);                                                               % 
Of=x(:,5);                                                              %
O=x(:,6);                                                               %
TR=x(:,7);                                                              %
TA=x(:,8);                                                              %
R=x(:,9);                                                                %
%
%
figure(1)
%
hold on;
plot(l,R,'r')
plot(l,C)
plot(l,Of,'g')
%plot(l,TR)
hold off;
xlabel('Time')
ylabel('Concetration (micromolar)')
title('Activator response')
%
end
%
%
function y=difeq(t,x)
%
global kcc kcd dsu duo fsu fuo gcs gcsa guf gu gof go gtr gta gr acs au ao ar atr ata btr bta br ho htr n
%
%--------------------------------------------%
%                               VARIABLES                                 %
%--------------------------------------------%
%
 C=input(t);                                    % Extracellular concentration of the cholerae autoinducer-1 (CAI-1)
%
 CS=x(1);                                       % Concentration of the membrane bound CqsS protein (CAI-1 unbound=Inactive)
 CSa=x(2);                                     % Concentration of the membrane bound CqsS protein (CAI-1 bound=Active)
 Uf= x(3);                                      % Phosphorelay protein LuxU (phosphorylated) Concentration
 U=x(4);                                        % Phosphorelay protein LuxU (unphosphorylated) Concentration
 Of=x(5);                                       % Phosphorelay protein LuxO (phosphorylated) Concentration
 O=x(6);                                        % Phosphorelay protein LuxO (unphosphorylated) Concentration
 TR=x(7);                                       % Ptet Repressor protein concentration
 TA=x(8);                                       % Ptet Activator protein concentration
 R=x(9);                                         % Response molecule concentration
% 
%--------------------------------------------%
%                               PARAMETERS                             %
%--------------------------------------------%
%
kcc=1;                                          % CAI1 and CqsS coupling Rate
kcd=1;                                          % CAI1 and CqsS decoupling Rate
dsu=4;                                          % LuxU* dephosphorylation rate through CqsS*
duo=4;                                         % LuxO* dephosphorylation rate through LuxU
fsu=2;                                          % LuxU* phosphorylation rate through CqsS
fuo=2;                                          % LuxO* phosphorylation rate through LuxU*
%
gcs=1;                                          % CqsS protein decay rate
gcsa=1;                                        % CqsS* protein decay rate
guf=1;                                          % LuxU* protein decay rate
gu=1;                                           % LuxU protein decay rate
gof=1;                                          % LuxO* protein decay rate
go=1;                                           % LuxO protein decay rate
gtr=1;                                          % Ptet Repressor protein decay rate
gta=1;                                          % Ptet Activator protein decay rate
gr=1;                                           % Response molecule decay rate
%
acs=3;                                          % CS basal production rate
au=3;                                           % LuxU basal production rate
ao=3;                                           % LuxO basal production rate
ar=0.01;                                       % response molecule basal production rate
atr=0.01;                                      % TR basal production rate
ata=0.01;                                      % TA basal production rate

btr=5;                                          % Maximum rate of TR expression
bta=5;                                          % Maximum rate of TA expression
br=5;                                           % Maximum rate of response molecule expression
ho=1.5;                                       % LuxO*- DNA coupling rate
htr=2;                                          % TRdomain-DNA coupling rate
%
n=1;                                            % Hill coefficient
%
%--------------------------------------------%
%                                 Equations                                %
%--------------------------------------------%
%
dCS = acs + kcd*CSa - C*CS*kcc - gcs*CS;                                      % Differential equation governing the change in the concentration of the membrane bound CqsS protein (CAI-1 unbound=Inactive)through time
dCSa = -kcd*CSa + C*CS*kcc - gcsa*CSa;                                        % Differential equation governing the change in the concentration of the membrane bound CqsS protein (CAI-1 bound=Active) through time
dUf = U*Of*duo-Uf*O*fuo-CSa*Uf*dsu + CS*U*fsu - guf*Uf;           % Differential equation governing the change in the phosphorelay protein LuxU (phosphorylated) Concentration through time
dU = au-U*Of*duo+Uf*O*fuo - CS*U*fsu + CSa*Uf*dsu -gu*U;       % Differential equation governing the change in the phosphorelay protein LuxU (unphosphorylated) Concentration through time
dOf = Uf*O*fuo - U*Of*duo - gof*Of;             				        % Differential equation governing the change in the phosphorelay protein LuxO (phosphorylated) Concentration through time
dO = ao - Uf*O*fuo + U*Of*duo - go*O;           		                % Differential equation governing the change in the phosphorelay protein LuxO (unphosphorylated) Concentration through time
dTR = atr + (btr*Of^n)/(ho^n+Of^n) - gtr*TR;   				% Differential equation governing the change in the ptet Repressor protein concentration through time
dTA = ata + bta/((1+TR/TA)+(htr/TA)) - gta*TA; 				% Differential equation governing the change in the ptet Activator protein concentration through time
dR = ar + br/((1+TR/TA)+(htr/TA)) - gr*R;       				% Differential equation governing the change in the response molecule concentration through time
%
y(1)=dCS;                                       %
y(2)=dCSa;                                     %
y(3)=dUf;                                       %
y(4)=dU;                                        %
y(5)=dOf;                                       %
y(6)=dO;                                        %
y(7)=dTR;                                       %
y(8)=dTA;                                       %
y(9)=dR;                                         %
%
y=y';                                              % Transposing the vector because of matlab language restrictions
%
end
%
%
function y=input(t)
%--------------------------------------------%
%                      CAI-1 (C) pulse simulation                 %
%--------------------------------------------%
if t>50 && t<100   
   y=10;
elseif t>250 && t<1500
    y=20;
else
    y=0;
end
%
%
</textarea>
<br>
<br>
<b><font color="#8A0808" size="3" >Metropolis-Hastings Algorithm</font></b>
<br>
This code determines which set of values for missing parameters yield the most desirable response.<br>
<textarea readonly rows = "20">
function y=model(patr, pbtr, pho, pkttr, pata, pbta, phtr,pn)
%
global kcc kcd kcu kuc kuo kou kttr gcs gcsa guf gu gof go gtr gta gr gttr acs au ao ar atr ata bcu buc btr buo bou bta ho htr n
%
% 
%--------------------------------------------%
%                              PARAMETERS                              %
%--------------------------------------------%
%
kcc=1;                                                  % CAI1 and CqsS coupling Rate
kcd=1;                                                  % CAI1 and CqsS decoupling Rate
kcu=2;                                                  % Michaelis Menten Constant for the phoshporylation CqsS-LuxU
kuc=1;                                                  % Michaelis Menten Constant of the phosphase process CqSS-LuxU
kuo=2;                                                  % Michaelis Menten Constant for the phoshporylation LuxU-LuxO
kou=1;                                                  % Michaelis Menten Constant of the phosphase process LuxU-LuxO
kttr=pkttr;                                             %Dimerization rate of TetR
%
gcs=0.5;                                                % CqsS protein decay rate
gcsa=1;                                                 % CqsS* protein decay rate
guf=0.12;                                              % LuxU* protein decay rate
gu=0.65;                                               % LuxU protein decay rate
gof=0.12;                                              % LuxO* protein decay rate
go=0.12;                                               % LuxO protein decay rate
gtr=0.12;                                              % Ptet Repressor protein decay rate
gttr=0.12;                                             %Ptet2 dimer Repressor decay rate
gta=0.12;                                              % Ptet Activator protein decay rate
gr=0.12;                                                % Response molecule decay rate
%
acs=5;                                                   % CS basal production rate
au=5;                                                    % LuxU basal production rate
ao=5;                                                    % LuxO basal production rate
ar=0;                                                     % response molecule basal production rate
atr=patr;                                                % TR basal production rate
ata=pata;                                               % TA basal production rate
%
bcu=1;                                                   % Phosphorylation rate of CqsS to LuxU
buc=0.1;                                                % Phosphase rate of CqsS to LuxU
buo=3.2;                                                % Phosphorylation rate LuxU-LuxO
bou=1.58;                                              % Phosphase rate LuxU-LuxO
btr=pbtr;                                                % Maximum rate of TR expression
bta=pbta;                                               % Maximum rate of TA expression
%
ho=pho;                                                 % LuxO*- DNA coupling rate
htr=phtr;                                                % TRdomain-DNA coupling rate
%
n=pn;                                                     % Hill coefficient
%
%
%
%
%
%--------------------------------------------%
%                Runge Kutta 4 order aproximation            %
%--------------------------------------------%
%
h=210;                                                                                  % Maximum time
m=1/12;                                                                               % Step lenght [h]
conInd=[1,0,1,1,1,1,0,1,1,1];                                                % 
l=(0:m:h)';                                                                             % Time vector
x=zeros(length(l),length(conInd));                                         % Variable's Matrix, the variables in each column and time in each row
C=zeros(1,length(l));                                                             % zero vector to be filled
x(1,:)=conInd;                                                                       % 
%
for k=1:length(l)-1
    xk=x(k,:);                                                                          % Capture of the last value in the matrix, the actual variable values
    k1=difeq_Metropolis(l(k),xk);                                             % First slope of the RK4 method
    k2=difeq_Metropolis(l(k)+m/2,xk+(m/2*k1)');                   % Second slope of the RK4 method
    k3=difeq_Metropolis(l(k)+m/2,xk+(m/2*k2)');                   % Third slope of the RK4 method
    k4=difeq_Metropolis(l(k)+m,xk+(m*k3)');                          % Fourth slope of the RK4 method
    xk1=xk+m/6*(k1+2*k2+2*k3+k4)';                                  % New value calculation for the variables
    %xk1=xk+m*difeq_Metropolis(l(k),xk)';                              % Another aproximation, Newton's method                      
    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
%
for j=1:length(l)
    if (l(j)>100 && l(j)<200)
        C(j)=20;
    else
        C(j)=0;
            
    end
end
%
    
CS=x(:,1);                                                             %
CSa=x(:,2);                                                           %
Uf=x(:,3);                                                             %
U=x(:,4);                                                              %
Of=x(:,5);                                                             %
O=x(:,6);                                                              % 
TR=x(:,7);                                                             %                      
TTR=x(:,8);                                                           %
TA=x(:,9);                                                             %
R=x(:,10);                                                             %

%figure;                                                                 %
%plot(l,C,l,R,'r');                                                     %
%legend ('CAI','Chromoprotein');                            %
%xlabel('Time [h]');                                                 %
%ylabel('Concentration [nanomolar]');                     %
%title('Response');                                                   %
%
y=[R(95*12),R(100*12),R(105*12),R(110*12),R(115*12),R(120*12),R(125*12),R(130*12),R(135*12),R(140*12),R(145*12),R(150*12),R(155*12),R(160*12),R(165*12),R(170*12)];
%
%
%
%
RFP=[3,3.1,3.15,3.5,4.8,6,11,15,16.5,18,18.5,19,19.5,20,20.25,20.5];
t=[95,100,105,110,115,120,125,130,135,140,145,150,155,160,165,170];
%
steps=5000;
%
atr_walk=zeros(1,steps+1);
btr_walk=zeros(1,steps+1);
ho_walk=zeros(1,steps+1);
kttr_walk=zeros(1,steps+1);
ata_walk=zeros(1,steps+1);
bta_walk=zeros(1,steps+1);
htr_walk=zeros(1,steps+1);
n_walk=zeros(1,steps+1);
%
atr_walk(1)=random('unif',0.5,1.5);
btr_walk(1)=random('unif',4,6);
ho_walk(1)=random('unif',1.2,1.8);
kttr_walk(1)=random('unif',2,7);
ata_walk(1)=random('unif',0.1,0.5);
bta_walk(1)=random('unif',4,7);
htr_walk(1)=random('unif',4,7);
n_walk(1)=random('unif',1,3);
%
    for i=1:steps
%
    atr_prime=random('unif',atr_walk(i)-0.2,atr_walk(i)+0.2);
    btr_prime=random('unif',btr_walk(i)-0.2,btr_walk(i)+0.2);
    ho_prime=random('unif',ho_walk(i)-0.2,ho_walk(i)+0.2);
    kttr_prime=random('unif',kttr_walk(i)-0.2,kttr_walk(i)+0.2);
    ata_prime=random('unif',ata_walk(i)-0.2,ata_walk(i)+0.2);
    bta_prime=random('unif',bta_walk(i)-0.2,bta_walk(i)+0.2);
    htr_prime=random('unif',htr_walk(i)-0.2,htr_walk(i)+0.2);
    n_prime=random('unif',n_walk(i)-0.2,n_walk(i)+0.2);
%  
    model_i=model(atr_walk(i),btr_walk(i),ho_walk(i),kttr_walk(i),ata_walk(i),bta_walk(i),htr_walk(i),n_walk(i));
    model_prime=model(atr_prime,btr_prime,ho_prime,kttr_prime,ata_prime,bta_prime,htr_prime,n_prime);
%
    alpha=likelihood_lineal(RFP,model_prime)/likelihood_lineal(RFP,model_i);
%    
        if(alpha>=1)
            atr_walk(i+1)=atr_prime;
            btr_walk(i+1)=btr_prime;
            ho_walk(i+1)=ho_prime;
            kttr_walk(i+1)=kttr_prime;
            ata_walk(i+1)=ata_prime;
            bta_walk(i+1)=bta_prime;
            htr_walk(i+1)=htr_prime;
            n_walk(i+1)=n_prime;
        else
           % beta=random('unif',0,1);
            %if(beta<=alpha)
             %  atr_walk(i+1)=atr_prime;
              % btr_walk(i+1)=btr_prime;
               %ho_walk(i+1)=ho_prime;
               %kttr_walk(i+1)=kttr_prime;
               %ata_walk(i+1)=ata_prime;
               %bta_walk(i+1)=bta_prime;
               %htr_walk(i+1)=htr_prime;
               %n_walk(i+1)=n_prime; 
            %else
               atr_walk(i+1)=atr_walk(i);
               btr_walk(i+1)=btr_walk(i);
               ho_walk(i+1)=ho_walk(i);
               kttr_walk(i+1)=kttr_walk(i);
               ata_walk(i+1)=ata_walk(i);
               bta_walk(i+1)=bta_walk(i);
               htr_walk(i+1)=htr_walk(i);
               n_walk(i+1)=n_walk(i); 
            %end
        end
    i
    end
%   
figure (1)
plot(t,RFP,'b',t,model(atr_walk(steps+1),btr_walk(steps+1),ho_walk(steps+1),kttr_walk(steps+1),ata_walk(steps+1),bta_walk(steps+1),htr_walk(steps+1),n_walk(steps+1)),'r')
xlabel('tiempo')
ylabel('Reportero')
title('Respuesta')
%
figure (2)
scatter((1:steps+1),atr_walk)
xlabel('step')
ylabel('valor')
title('atr')
%
figure (3)
scatter((1:steps+1),btr_walk)
xlabel('step')
ylabel('valor')
title('btr')
%
figure (4)
scatter((1:steps+1),ho_walk)
xlabel('step')
ylabel('valor')
title('ho')
%
figure (5)
scatter((1:steps+1),kttr_walk)
xlabel('step')
ylabel('valor')
title('kttr')
%
figure (6)
scatter((1:steps+1),ata_walk)
xlabel('step')
ylabel('valor')
title('ata')
%
figure (7)
scatter((1:steps+1),bta_walk)
xlabel('step')
ylabel('valor')
title('bta')
%
figure (8)
scatter((1:steps+1),htr_walk)
xlabel('step')
ylabel('valor')
title('htr')
%
figure (9)
scatter((1:steps+1),n_walk)
xlabel('step')
ylabel('valor')
title('n')
%
%
function y=input_Metropolis(t)
%--------------------------------------------%
%                    CAI-1 (C) pulse simulation                    %
%--------------------------------------------%
if t>100 && t<200   
   y=20;
else
    y=0;
end
</textarea>
<br>
<br>
<b><font color="#8A0808" size="3" >Sensitivity Analysis</font></b>
<br>
This code points out which are the critical parameters in the system (those that change the response drastically).<br>
<textarea readonly rows = "20">
%
global kcc kcd kcu kuc kuo kou kttr gcs gcsa guf gu gof go gtr gta gr gttr acs au ao ar atr ata bcu buc btr buo bou bta ho htr n
%
% 
%--------------------------------------------%
%                                PARAMETERS                            %
%--------------------------------------------%
%
kcc=1;                                                           % CAI1 and CqsS coupling Rate
kcd=1;                                                           % CAI1 and CqsS decoupling Rate
kcu=2;                                                           % Michaelis Menten Constant for the phoshporylation CqsS-LuxU
kuc=1;                                                           % Michaelis Menten Constant of the phosphase process CqSS-LuxU
kuo=2;                                                           % Michaelis Menten Constant for the phoshporylation LuxU-LuxO
kou=1;                                                           % Michaelis Menten Constant of the phosphase process LuxU-LuxO
kttr=5;                                                           % Dimerization rate of TetR ---[Sensitivity Analysis]--- [5-50/5]
%
gcs=0.5;                                                         % CqsS protein decay rate
gcsa=1;                                                          % CqsS* protein decay rate
guf=0.12;                                                       % LuxU* protein decay rate
gu=0.65;                                                        % LuxU protein decay rate
gof=0.12;                                                       % LuxO* protein decay rate
go=0.12;                                                        % LuxO protein decay rate
gtr=0.12;                                                       % Ptet Repressor protein decay rate
gttr=0.12;                                                      % Ptet2 dimer Repressor decay rate
gta=0.12;                                                       % Ptet Activator protein decay rate
gr=0.12;                                                         % Response molecule decay rate
%
acs=5;                                                             % CS basal production rate
au=5;                                                              % LuxU basal production rate
ao=5;                                                              % LuxO basal production rate
ar=0;                                                               % response molecule basal production rate
atr=1;                                                              % TR basal production rate  ---[Sensitivity Analysis]--- [0.5-5/0.5]
ata=0.1;                                                           % TA basal production rate ---[Sensitivity Analysis]--- [0.1-10/0.5]
%
bcu=1;                                                             % Phosphorylation rate of CqsS to LuxU
buc=0.1;                                                          % Phosphase rate of CqsS to LuxU
buo=3.2;                                                          % Phosphorylation rate LuxU-LuxO
bou=1.58;                                                        % Phosphase rate LuxU-LuxO
btr=5;                                                              % Maximum rate of TR expression ---[Sensitivity Analysis]--- [1-10/1]
bta=5;                                                              % Maximum rate of TA expression ---[Sensitivity Analysis]--- [1-10/1]
%
ho=1.5;                                                            % LuxO*- DNA coupling rate ---[Sensitivity Analysis]--- [0.5-2.5/0.5]
htr=5;                                                               % TRdomain-DNA coupling rate ---[Sensitivity Analysis]--- [1-10/1]
%
n=2;                                                                  % Hill's coefficient ---[Sensitivity Analysis]--- [1-4/0.5]
%
h=400;                                                              % Maximum time
m=1/10;                                                            % Step lenght [h]
%
param= (5:50:5);                                                % Range for the parameter to vary, evaluate at least 2 orders of magnitude chose maximum 20 points.
rta=zeros(3,length(param));                                % Array that saves the canche in concentration of LuxOF, the response and time of response.
for w=1:length(param)                                        % Varying the value of the parameter
    kttr= param(w);                                               % New parameter value (This MUST be changed everytime a parameter is being evaluated)
%
%
%--------------------------------------------%
%               Runge Kutta 4 order aproximation             %
%--------------------------------------------%
%
%
%
conInd=[1,0,1,1,1,1,0,1,1,1];                               %
l=(0:m:h)';                                                            % Time vector
%
x=zeros(length(l),length(conInd));                         % Variable's Matrix, the variables in each column and time in each row
C=zeros(1,length(l));                                             % zero vector to be filled
x(1,:)=conInd;                                                       %
%
for k=1:length(l)-1
    xk=x(k,:);                                                           % Capture of the last value in the matrix, the actual variable values
    k1=difeq(l(k),xk);                                               % First slope of the RK4 method
    k2=difeq(l(k)+m/2,xk+(m/2*k1)');                     % Second slope of the RK4 method
    k3=difeq(l(k)+m/2,xk+(m/2*k2)');                     % Third slope of the RK4 method
    k4=difeq(l(k)+m,xk+(m*k3)');                            % Fourth slope of the RK4 method
    xk1=xk+m/6*(k1+2*k2+2*k3+k4)';                  % New value calculation 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(k+1,:)=xk2;                                % Actualizacion del nuevo vector de variables en la matriz
end
%
 R=x(:,10); 
 Of=x(:,5);%
%
%
 it=141;
%
  for it=1011:3001 
%
%
   if (R(it)-R(1011))/(R(1011))> 1.5      % If the response changes 1.5-fold relative to previous iteration
%
   rta(1,w)= l(it)-100;
%
%
%
  end
%
%
%
  end   
%
%
%
%
rta(2,w)= Of(2701)-Of(801);
rta(3,w)= R(2701)-R(801);
%
end
%
%Here Labels should be changed for them to make proper sense for each evaluated parameter
%
figure (1)                               
plot(param,rta(1,:),'r')
xlabel('Dimerization rate of TetR');                                    
ylabel('Time of response')
title('Sensitivity analysis')
%
figure (2)
plot(param,rta(2,:),'g')
xlabel('Dimerization rate of TetR');
ylabel('Change in LuxO Phosphorylated')
title('Sensitivity analysis')
%
figure (3)
plot(param,rta(3,:))   
xlabel('Dimerization rate of TetR');
ylabel('Change in response (nM)')
title('Sensitivity analysis')
</textarea>
<br>
<br>
<b><font color="#8A0808" size="3" >Stochastic Model</font></b>
<br>
Sometimes probabilistic models better describe certain systems, This required a fairly big amount of computational power<br>
<textarea readonly rows = "20">
%Stochastic Simulation for the system
%
%%----Parameters of the System----
global kcc kcd kcu kuc kuo kou kttr gcs gcsa guf gu gof go gtr gta gr acs au ao ar atr ata bcu buc btr buo bou bta ho htr n gttr
%
%
numCells = 20;
numEvents= 50000;
%
%
%
% 
kcc=1;                                               % CAI1 and CqsS coupling Rate
kcd=1;                                                % CAI1 and CqsS decoupling Rate
kcu=1.2044000;                                        % Michaelis Menten Constant for the phoshporylation CqsS-LuxU
kuc=0.6022000;                                          % Michaelis Menten Constant of the phosphase process CqSS-LuxU
kuo=1.2044000;                                         % Michaelis Menten Constant for the phoshporylation LuxU-LuxO
kou=0.6022000;                                        % Michaelis Menten Constant of the phosphase process LuxU-LuxO
kttr=5;                                              %Dimerization rate of TetR
%
gcs=0.500;                                              % CqsS protein decay rate
gcsa=1;                                             % CqsS* protein decay rate
guf=0.120;                                              % LuxU* protein decay rate
gu=0.650;                                               % LuxU protein decay rate
gof=0.120;                                              % LuxO* protein decay rate
go=0.650;                                               % LuxO protein decay rate
gtr=0.120;                                              % Ptet Repressor protein decay rate
gttr=0.120;                                              %Ptet2 dimer Repressor decay rate
gta=0.120;                                              % Ptet Activator protein decay rate
gr=0.050;                                               % Response molecule decay rate
%
acs=5;                                              % CS basal production rate
au=5;                                               % LuxU basal production rate
ao=5;                                               % LuxO basal production rate
ar=0.1;                                            % response molecule basal production rate
atr=1;                                           % TR basal production rate
ata=0.1;                                           % TA basal production rate
%
bcu=1;                                              % Phosphorylation rate of CqsS to LuxU
buc=0.1;                                              % Phosphase rate of CqsS to LuxU
buo=3.2;                                              % Phosphorylation rate LuxU-LuxO
bou=1.58;                                              % Phosphase rate LuxU-LuxO
btr=5;                                              % Maximum rate of TR expression
bta=5;                                              % Maximum rate of TA expression
%
ho=1.5;                                               % LuxO*- DNA coupling rate
htr=5;                                              % TRdomain-DNA coupling rate
%
n=3;% Hill coefficient
%
%%----Stochastic Simulation----
%
%Variables of the System
C=zeros(numCells,numEvents);
CS=zeros(numCells,numEvents);
CSa=zeros(numCells,numEvents);
Uf=zeros(numCells,numEvents);
U=zeros(numCells,numEvents);
Of=zeros(numCells,numEvents);
O=zeros(numCells,numEvents);
TR=zeros(numCells,numEvents);
TRR=zeros(numCells,numEvents);
TA=zeros(numCells,numEvents);
R=zeros(numCells,numEvents);
%
C(:,1)=0;
CS(:,1)=10;
CSa(:,1)=1;
Uf(:,1)=15;
U(:,1)=5;
Of(:,1)=40;
O(:,1)=5;
TA(:,1)=10;
TR(:,1)=2;
TRR(:,1)=1;
%
%
%
%
t = zeros(numCells,numEvents); %Time
tr = zeros(numCells,numEvents);
%
%
%Stochastic Simulation
%
for i =1:numCells
    
    i
    for j = 1:numEvents
        
        j
        
        %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= kcc*CS(i,j)*C(i,j);
        e2= kcd*CSa(i,j);
        e3= gcs*CS(i,j);
        e4= acs;
        e5= gcsa*CSa(i,j);
        e6= au;
        e7= (bcu*U(i,j)*CS(i,j))/(kcu+U(i,j));
        e8= (buc*Uf(i,j)*CSa(i,j))/(kuc+Uf(i,j));
        e9= (buo*O(i,j)*Uf(i,j))/(kuo+O(i,j));
        e10=(bou*Of(i,j)*U(i,j))/(kou+Of(i,j));
        e11=gu*U(i,j);
        e12=guf*Uf(i,j);
        e13=ao;
        e14=go*O(i,j);
        e15=gof*Of(i,j);
        
        e16=atr;
        e17=(btr*Of(i,j)^n)/((ho.^n)+(Of(i,j).^n));
        e18=kttr*TR(i,j)^2;
        e19=gtr*TR(i,j);
        
        e20=ata;
        e21=ar;
        e22=bta/(1+(TRR(i,j)/TA(i,j))+(htr/TA(i,j)));
        e23=gta*TA(i,j);
        e24=gr*R(i,j);
        e25=TRR(i,j)*gttr;
        
        
        %Vector of Events
        vecEvents = [e1,e2,e3,e4,e5,e6,e7,e8,e9,e10,e11,e12,e13,e14,e15,e16,e17,e18,e19,e20,e21,e22,e23,e24,e25];
        
        %Sum of Events
        sumEvents = sum(vecEvents);
        
        %Waiting time until the next ocurrence of a reaction time
        %Probabilistic Distribution
        tau= (1/sumEvents)*log(1/rand1);
        tr(i,j)=tau;
        t(i,j+1)= t(i,j)+ tau;
        
        %CAI outside the cell
        
        C(i,j+1)=funcImpulso(t(i,j+1));
        
        %if(t(i,j)<130)
         %   C(i,j)=funcImpulso(t(i,j+1))*957/2575;
        %end
        
        %Simulation Creation/Destruction of each molecule
        
        if (rand2<(vecEvents(1)/sumEvents))  %CqsS activation
            
            if (CS(i,j)<=0)
                
                CS(i,j+1)=0;
            else
                CS(i,j+1)= CS(i,j)-1;
            end
            
            CSa(i,j+1)= CSa(i,j)+1;
            U(i,j+1)= U(i,j);
            Uf(i,j+1)= Uf(i,j);
            O(i,j+1)= O(i,j);
            Of(i,j+1)= Of(i,j);
            TR(i,j+1)= TR(i,j);
            TRR(i,j+1)= TRR(i,j);
            TA(i,j+1)= TA(i,j);
            R(i,j+1)= R(i,j);
            
            
            
        elseif (rand2< (sum(vecEvents(1:2)))/sumEvents) %CqsS desactivation
            
            if (CSa(i,j)<=0)
                
                CSa(i,j+1)= 0;
                
            else
                CSa(i,j+1)= CSa(i,j)-1;
            end
            
            CS(i,j+1)= CS(i,j)+1;
            U(i,j+1)= U(i,j);
            Uf(i,j+1)= Uf(i,j);
            O(i,j+1)= O(i,j);
            Of(i,j+1)= Of(i,j);
            TR(i,j+1)= TR(i,j);
            TRR(i,j+1)= TRR(i,j);
            TA(i,j+1)= TA(i,j);
            R(i,j+1)= R(i,j);
            
            
            
        elseif (rand2< (sum(vecEvents(1:3)))/sumEvents) %CqsS degradation
            
            if (CS(i,j)<=0)
                
                CS(i,j+1)=0;
            else
                
                CS(i,j+1)= CS(i,j)-1;
            end
            
            CSa(i,j+1)= CSa(i,j);
            U(i,j+1)= U(i,j);
            Uf(i,j+1)= Uf(i,j);
            O(i,j+1)= O(i,j);
            Of(i,j+1)= Of(i,j);
            TR(i,j+1)= TR(i,j);
            TRR(i,j+1)= TRR(i,j);
            TA(i,j+1)= TA(i,j);
            R(i,j+1)= R(i,j);
            
        elseif (rand2< (sum(vecEvents(1:4)))/sumEvents)  %CqsS basal production
            
            
            CS(i,j+1)= CS(i,j)+1;
            CSa(i,j+1)= CSa(i,j);
            U(i,j+1)= U(i,j);
            Uf(i,j+1)= Uf(i,j);
            O(i,j+1)= O(i,j);
            Of(i,j+1)= Of(i,j);
            TR(i,j+1)= TR(i,j);
            TRR(i,j+1)= TRR(i,j);
            TA(i,j+1)= TA(i,j);
            R(i,j+1)= R(i,j);
            
        elseif (rand2< (sum(vecEvents(1:5)))/sumEvents) %CqsS* degradation
            
            if (CSa(i,j)<=0)
                
                CSa(i,j+1)=0;
            else
                
                CSa(i,j+1)= CSa(i,j)-1;
            end
            
            CS(i,j+1)= CS(i,j);
            U(i,j+1)= U(i,j);
            Uf(i,j+1)= Uf(i,j);
            O(i,j+1)= O(i,j);
            Of(i,j+1)= Of(i,j);
            TR(i,j+1)= TR(i,j);
            TRR(i,j+1)= TRR(i,j);
            TA(i,j+1)= TA(i,j);
            R(i,j+1)= R(i,j);
            
            
        elseif (rand2< (sum(vecEvents(1:6)))/sumEvents)  %LuxU basal production
            
            
            U(i,j+1)= U(i,j)+1;
            CS(i,j+1)= CS(i,j);
            CSa(i,j+1)= CSa(i,j);
            Uf(i,j+1)= Uf(i,j);
            O(i,j+1)= O(i,j);
            Of(i,j+1)= Of(i,j);
            TR(i,j+1)= TR(i,j);
            TRR(i,j+1)= TRR(i,j);
            TA(i,j+1)= TA(i,j);
            R(i,j+1)= R(i,j);
            
            
        elseif (rand2< (sum(vecEvents(1:7)))/sumEvents)  %LuxU fosforilation
            
            if (U(i,j)<=0)
                
                U(i,j+1)=0;
            else
                U(i, j+1)= U(i,j)-1;
            end
            
            CS(i,j+1)= CS(i,j);
            CSa(i,j+1)= CSa(i,j);
            Uf(i,j+1)= Uf(i,j)+1;
            O(i,j+1)= O(i,j);
            Of(i,j+1)= Of(i,j);
            TR(i,j+1)= TR(i,j);
            TRR(i,j+1)= TRR(i,j);
            TA(i,j+1)= TA(i,j);
            R(i,j+1)= R(i,j);
            
            
        elseif (rand2< (sum(vecEvents(1:8)))/sumEvents)  %LuxU dephosphorilation
            
            if (Uf(i,j)<=0)
                
                Uf(i,j+1)=0;
            else
                Uf(i, j+1)= Uf(i,j)-1;
            end
            
            CS(i,j+1)= CS(i,j);
            CSa(i,j+1)= CSa(i,j);
            U(i,j+1)= U(i,j)+1;
            O(i,j+1)= O(i,j);
            Of(i,j+1)= Of(i,j);
            TR(i,j+1)= TR(i,j);
            TRR(i,j+1)= TRR(i,j);
            TA(i,j+1)= TA(i,j);
            R(i,j+1)= R(i,j);
            
            
        elseif (rand2< (sum(vecEvents(1:9)))/sumEvents)  %LuxO phosphorilation
            
            if (Uf(i,j)<=0 || O(i,j)<=0 )
                if(Uf(i,j)<0)
                    Uf(i,j+1)=0;
                else
                    Uf(i,j+1)=Uf(i,j);
                end
                if(O(i,j)<0)
                    O(i,j+1)=0;
                else
                    O(i,j+1)=O(i,j);
                end
              
            else
                Uf(i, j+1)= Uf(i,j)-1;
                O(i,j+1)=O(i,j)-1;
            end
            
            CS(i,j+1)= CS(i,j);
            CSa(i,j+1)= CSa(i,j);
            U(i,j+1)= U(i,j)+1;
            Of(i,j+1)= Of(i,j)+1;
            TR(i,j+1)= TR(i,j);
            TRR(i,j+1)= TRR(i,j);
            TA(i,j+1)= TA(i,j);
            R(i,j+1)= R(i,j);
            
        elseif (rand2< (sum(vecEvents(1:10)))/sumEvents)  %LuxO dephosphorilation
            
            if (U(i,j)<=0)
                
                U(i,j+1)=0;
            else
                U(i,j+1)= U(i,j)-1;
            end
            
            if (Of(i,j)<=0)
                
                Of(i,j+1)=0;
            else
                Of(i, j+1)= Of(i,j)-1;
            end
            
            
            CS(i,j+1)= CS(i,j);
            CSa(i,j+1)= CSa(i,j);
            Uf(i,j+1)= Uf(i,j)+1;
            O(i,j+1)= O(i,j)+1;
            TR(i,j+1)= TR(i,j);
            TRR(i,j+1)= TRR(i,j);
            TA(i,j+1)= TA(i,j);
            R(i,j+1)= R(i,j);
            
        elseif (rand2< (sum(vecEvents(1:11)))/sumEvents)  %LuxU degradation
            
            if (U(i,j)<=0)
                
                U(i,j+1)=0;
            else
                U(i, j+1)= U(i,j)-1;
            end
            
            CS(i,j+1)= CS(i,j);
            CSa(i,j+1)= CSa(i,j);
            Uf(i,j+1)= Uf(i,j);
            O(i,j+1)= O(i,j);
            Of(i,j+1)= Of(i,j);
            TR(i,j+1)= TR(i,j);
            TRR(i,j+1)= TRR(i,j);
            TA(i,j+1)= TA(i,j);
            R(i,j+1)= R(i,j);
       
        elseif (rand2< (sum(vecEvents(1:12)))/sumEvents )  %LuxU* degradation 
            
            if (Uf(i,j)<=0)
                
                Uf(i,j+1)=0;
            else
                Uf(i, j+1)= Uf(i,j)-1;
            end
            
            CS(i,j+1)= CS(i,j);
            CSa(i,j+1)= CSa(i,j);
            U(i,j+1)= U(i,j);
            O(i,j+1)= O(i,j);
            Of(i,j+1)= Of(i,j);
            TR(i,j+1)= TR(i,j);
            TRR(i,j+1)= TRR(i,j);
            TA(i,j+1)= TA(i,j);
            R(i,j+1)= R(i,j);
            
            
        elseif (rand2< (sum(vecEvents(1:13)))/sumEvents )%Basal production LuxO 
            
            
            U(i, j+1)= U(i,j);
            CS(i,j+1)= CS(i,j);
            CSa(i,j+1)= CSa(i,j);
            Uf(i,j+1)= Uf(i,j);
            O(i,j+1)= O(i,j)+1;
            Of(i,j+1)= Of(i,j);
            TR(i,j+1)= TR(i,j);
            TRR(i,j+1)= TRR(i,j);
            TA(i,j+1)= TA(i,j);
            R(i,j+1)= R(i,j);
            
            
        elseif (rand2< (sum(vecEvents(1:14)))/sumEvents) %LuxO degradation 
            
            if (O(i,j)<=0)
                
                O(i,j+1)=0;
            else
                O(i, j+1)= O(i,j)-1;
            end
            
            CS(i,j+1)= CS(i,j);
            CSa(i,j+1)= CSa(i,j);
            U(i,j+1)= U(i,j);
            Uf(i,j+1)= Uf(i,j);
            Of(i,j+1)= Of(i,j);
            TR(i,j+1)= TR(i,j);
            TRR(i,j+1)= TRR(i,j);
            TA(i,j+1)= TA(i,j);
            R(i,j+1)= R(i,j);
            
        elseif (rand2< (sum(vecEvents(1:15)))/sumEvents ) %LuxO* degradation 
            
            if (Of(i,j)<=0)
                
                Of(i,j+1)=0;
            else
                Of(i, j+1)= Of(i,j)-1;
            end
            
            CS(i,j+1)= CS(i,j);
            CSa(i,j+1)= CSa(i,j);
            U(i,j+1)= U(i,j);
            O(i,j+1)= O(i,j);
            Uf(i,j+1)= Uf(i,j);
            TR(i,j+1)= TR(i,j);
            TRR(i,j+1)= TRR(i,j);
            TA(i,j+1)= TA(i,j);
            R(i,j+1)= R(i,j);
            
        elseif (rand2< (sum(vecEvents(1:16)))/sumEvents) %Basal Production tetR 
            
            
            TR(i, j+1)= TR(i,j)+1;
            CS(i,j+1)= CS(i,j);
            CSa(i,j+1)= CSa(i,j);
            U(i,j+1)= U(i,j);
            O(i,j+1)= O(i,j);
            Of(i,j+1)= Of(i,j);
            Uf(i,j+1)= Uf(i,j);
            TRR(i,j+1)= TRR(i,j);
            TA(i,j+1)= TA(i,j);
            R(i,j+1)= R(i,j);
            
        elseif (rand2< (sum(vecEvents(1:17)))/sumEvents )%Active production tetR 
            
            
            TR(i, j+1)= TR(i,j)+1;
            CS(i,j+1)= CS(i,j);
            CSa(i,j+1)= CSa(i,j);
            U(i,j+1)= U(i,j);
            O(i,j+1)= O(i,j);
            Of(i,j+1)= Of(i,j);
            Uf(i,j+1)= Uf(i,j);
            TRR(i,j+1)= TRR(i,j);
            TA(i,j+1)= TA(i,j);
            R(i,j+1)= R(i,j);
           
            
        elseif (rand2< (sum(vecEvents(1:18)))/sumEvents) %tetR dimerization  
            
            if (TR(i,j)<=0)
                
                TR(i,j+1)=0;
                TRR(i, j+1)= TRR(i,j);
                
            elseif (TR(i,j)==1)
                TR(i, j+1)=1;
                TRR(i, j+1)= TRR(i,j);
            else
                TR(i, j+1)= TR(i,j)-2;
                TRR(i, j+1)= TRR(i,j)+1;
            end
            
            CS(i,j+1)= CS(i,j);
            CSa(i,j+1)= CSa(i,j);
            U(i,j+1)= U(i,j);
            O(i,j+1)= O(i,j);
            Of(i,j+1)= Of(i,j);
            Uf(i,j+1)= Uf(i,j);
            TA(i,j+1)= TA(i,j);
            R(i,j+1)= R(i,j);
            
        elseif (rand2< (sum(vecEvents(1:19)))/sumEvents) %tetR degradation  
            
            if (TR(i,j)<=0)
                
                TR(i,j+1)=0;
            else
                TR(i,j+1)= TR(i,j)-1;
            end
            
            CS(i,j+1)= CS(i,j);
            CSa(i,j+1)= CSa(i,j);
            U(i,j+1)= U(i,j);
            O(i,j+1)= O(i,j);
            Of(i,j+1)= Of(i,j);
            Uf(i,j+1)= Uf(i,j);
            TRR(i,j+1)= TRR(i,j);
            TA(i,j+1)= TA(i,j);
            R(i,j+1)= R(i,j);
        
        elseif (rand2< (sum(vecEvents(1:20)))/sumEvents )%Basal production Tet activator  
            
            
            TA(i, j+1)= TA(i,j)+1;
            CS(i,j+1)= CS(i,j);
            CSa(i,j+1)= CSa(i,j);
            U(i,j+1)= U(i,j);
            O(i,j+1)= O(i,j);
            Of(i,j+1)= Of(i,j);
            Uf(i,j+1)= Uf(i,j);
            TRR(i,j+1)= TRR(i,j);
            TR(i,j+1)= TR(i,j);
            R(i,j+1)= R(i,j);
        
        elseif (rand2< (sum(vecEvents(1:21)))/sumEvents) %Basal production Reporter protein  
            
  
            R(i, j+1)= R(i,j)+1;
            CS(i,j+1)= CS(i,j);
            CSa(i,j+1)= CSa(i,j);
            U(i,j+1)= U(i,j);
            O(i,j+1)= O(i,j);
            Of(i,j+1)= Of(i,j);
            Uf(i,j+1)= Uf(i,j);
            TRR(i,j+1)= TRR(i,j);
            TR(i,j+1)= TR(i,j);
            TA(i,j+1)= TA(i,j);
            
        elseif (rand2< (sum(vecEvents(1:22)))/sumEvents )%Active production Reporter protein & activator  
            
            TA(i,j+1)= TA(i,j)+1;
            R(i,j+1)= R(i,j)+1;
            CS(i,j+1)= CS(i,j);
            CSa(i,j+1)= CSa(i,j);
            U(i,j+1)= U(i,j);
            O(i,j+1)= O(i,j);
            Of(i,j+1)= Of(i,j);
            Uf(i,j+1)= Uf(i,j);
            TRR(i,j+1)= TRR(i,j);
            TR(i,j+1)= TR(i,j);
            
        elseif (rand2< (sum(vecEvents(1:23)))/sumEvents) %Tet activator degradation  
            
            if (TA(i,j)<=0)
                
                TA(i,j+1)=0;
            else
                TA(i, j+1)= TA(i,j)-1;
            end
            
            CS(i,j+1)= CS(i,j);
            CSa(i,j+1)= CSa(i,j);
            U(i,j+1)= U(i,j);
            O(i,j+1)= O(i,j);
            Of(i,j+1)= Of(i,j);
            Uf(i,j+1)= Uf(i,j);
            TRR(i,j+1)= TRR(i,j);
            TR(i,j+1)= TR(i,j);
            R(i,j+1)= R(i,j);
            
        elseif (rand2 < (sum(vecEvents(1:24)))/sumEvents) %Reporter protein activator degradation  
            
            if (R(i,j)<=0)
                R(i,j+1)=0;
            else
                R(i, j+1)= R(i,j)-1;
            end
            
            CS(i,j+1)= CS(i,j);
            CSa(i,j+1)= CSa(i,j);
            U(i,j+1)= U(i,j);
            O(i,j+1)= O(i,j);
            Of(i,j+1)= Of(i,j);
            Uf(i,j+1)= Uf(i,j);
            TRR(i,j+1)= TRR(i,j);
            TR(i,j+1)= TR(i,j);
            TA(i,j+1)= TA(i,j);
            
        
        elseif (rand2 < (sum(vecEvents(1:25)))/sumEvents) %TetR dimers degradation  
            
            if (TRR(i,j)<=0)
                TRR(i,j+1)=0;
            else
                TRR(i, j+1)= TRR(i,j)-1;
            end
            
            CS(i,j+1)= CS(i,j);
            CSa(i,j+1)= CSa(i,j);
            U(i,j+1)= U(i,j);
            O(i,j+1)= O(i,j);
            Of(i,j+1)= Of(i,j);
            Uf(i,j+1)= Uf(i,j);
            TR(i,j+1)= TR(i,j);
            TA(i,j+1)= TA(i,j);
            R(i, j+1)= R(i,j);
            
        end    
    end
end
%
Pro_CS=zeros(numEvents);
Pro_CSa=zeros(numEvents);
Pro_Uf=zeros(numEvents);
Pro_U=zeros(numEvents);
Pro_Of=zeros(numEvents);
Pro_O=zeros(numEvents);
Pro_TR=zeros(numEvents);
Pro_TRR=zeros(numEvents);
Pro_TA=zeros(numEvents);
Pro_R=zeros(numEvents);
%
%
for i=1:numEvents
   Pro_CS(i)=(sum(CS(:,i)))/numCells;
end
for i=1:numEvents
    Pro_CSa(i)=(sum(CSa(:,i)))/numCells;
end
for i=1:numEvents
    Pro_Uf(i)=(sum(Uf(:,i)))/numCells;
end
for i=1:numEvents
    Pro_U(i)=(sum(U(:,i)))/numCells;
end
for i=1:numEvents
   Pro_Of(i)=(sum(Of(:,i)))/numCells;
end
for i=1:numEvents
   Pro_O(i)=(sum(O(:,i)))/numCells;
end
for i=1:numEvents
   Pro_TR(i)=(sum(TR(:,i)))/numCells;
end
for i=1:numEvents
   Pro_TRR(i)=(sum(TRR(:,i)))/numCells;
end
for i=1:numEvents
   Pro_TA(i)=(sum(TA(:,i)))/numCells;
end
for i=1:numEvents
   Pro_R(i)=(sum(R(:,i)))/numCells;
end
%
%
%
function [newTime newData] = regIntervalFixed2(Time,Data,step,timeMax)
%
% This function takes the output of a gillespie's simulation and turns it
% into a matrix with a time vector of regular intervals.
%
% INPUT
% Time: A  matrix sized number of Simulations(Cells) by number of Events
%   that contains the time asociated with each event in each simulation.
% Data: A matrix with the data. (events in columns and simulation in rows)
% step: a scalar of the time interval. It must have the same units as the
%   entrances of the matrix Time.
% timeMax: the maximum time required. It must have the same units as the
%   entrances of the matrix Time.
%
% OUTPUT:
% newTime: a vector of time regularly spaced. 
% newData: a matrix where each row represents a simulation (cell) and each
%   column represents the number of molecules for each time in newTime

L=size(Time,2);

newTime=0:step:timeMax;
newData = zeros(size(Data,1),length(newTime));

for k=1:size(Data,1)
    newdata=zeros(size(newTime));  
    j=1;
    for i=1:length(newTime)

      if step*(i-1)>Time(k,j) && j<L
          while step*(i-1)>Time(k,j) && j<L
              j=j+1;
          end    
      end
      newdata(i)=Data(k,j);
    end
    newData(k,:) = newdata;
end
%newdata(L+1)=data(L);
%
%
function impulse= funcImpulso(t)

if (t>150 && t<310)
    impulse=100; % Molecules/L
else
    impulse=0; %Molecules/L
end
%
end
%
%
function w=condiciones()
%
xi=[0.5 1000 5 1000];
%
y=fsolve(@CondInd,xi,optimset('algorithm','levenberg-marquardt','maxiter',100000,'tolfun',1e-9));
%
w=y;
%
end
%
%
</textarea>
</p>
</div>
<br>
<br>
<br>
<br>
<br>


<center>
<div class="button-fill orange" ><div class="button-text"> Back to modeling</div><div class="button-inside"><div class="inside-text"><a style="text-decoration: none; background-color: none; color: red;" href="https://2014.igem.org/Team:Colombia/Modeling">Go! </a></div></div></div>
</center>
<script>
     $(".button-fill").hover(function () {
    $(this).children(".button-inside").addClass('full');
}, function() {
  $(this).children(".button-inside").removeClass('full');
});
    //@ sourceURL=pen.js
  </script>

</html>