Team:Colombia/Modeling/Scripting

From 2012.igem.org

(Difference between revisions)
(→OPTIMIZATION)
(→STOCHASTIC MODEL)
Line 924: Line 924:
==STOCHASTIC MODEL ==
==STOCHASTIC MODEL ==
 +
 +
% MATLAB Code for modeling a gene net using Gillespie algorithm
 +
% Code by Roberto Moran, Daniela Olivera, Cesar Quintana & Andrés Simbaqueba
 +
% iGEM Colombia Team 2012, Universidad de Los Andes, Bogotá (Colombia)
 +
%
 +
%-------------------------------------------------------------------------%
 +
%                                  Parameters
 +
%-------------------------------------------------------------------------%
 +
clear all;
 +
%
 +
conv=602.2; %Conversion factor from micromolar to molecules/cell
 +
%
 +
%
 +
alfA = 0.9*conv;      %Basal production of Chitinase inside the cell (micromolar)
 +
alfP = 0.9*conv;      %Basal production of chitoporin
 +
alfC = 0.9*conv;      %Basal production of the CBP
 +
alfR = 0.06*conv;      %Basal production of LuxR
 +
alfI = 0.04*conv;      %Basal production of LuxI
 +
alfCI = 0.01*conv;    %Basal production of CI
 +
alfHA = 0.7*conv;    %Basal production of HipA7
 +
alfHB = 0.2*conv;    %Basal production of HipB
 +
alfAS = 0.4*conv;    %Basal production of Salycilic acid
 +
%
 +
gammaA=1;    %Degradation of Chitinase inside the cell
 +
gammaP=1;    %Degradation of chitoporin
 +
gammaC=1;    %Degradation concentration of the CBP
 +
gammaR=1;    %Degradation of LuxR
 +
gammaI=1;    %Degradation of LuxI
 +
gammaCI=4;  %Degradation of CI
 +
gammaHA=1;  %Degradation of HipA7
 +
gammaHB=4;  %Degradation of HipB
 +
gammaAS=0.8;  %Degradation of Salycilic acid
 +
gammaCS=1;  %Degradation of the complex Cs
 +
%
 +
mCS=13/conv^2;      %Kinetic constant for the formation of the complex CS
 +
mCSQ=12/conv^2;    %Kinetic constant of the reaction of the complex CS with the chitin
 +
mAQQ=0.2/conv^2;    %Kinetic constant for the reaction of the chitinase and th chitin
 +
mIR=3/conv^2;      %Kinetic constant for the formation of the complex LuxILuxR
 +
mI=3/conv;      %Constant that represent the union of the complex LuxILuxR with the promoter
 +
mHAHB=12/conv^4;    %Kinetic constant for the inhibition of HipA7
 +
%
 +
betaP=12*conv;    %Max production of the chitoporin
 +
betaA=10*conv;    %Max production of chitinase
 +
betaI=1*conv;    %Max production of LuxI
 +
betaCI=1*conv;%9.96*conv;  %Max production of CI
 +
betaHB=9.5*conv;  %Max production of HipB
 +
betaHA=11*conv;  %Max production of HipA7
 +
betaAS=1.12*conv;%11.2*conv;  %Max production of Salicylic acid
 +
%
 +
kS=0.08*conv;      %Constant k of the hill ecuation for the promoter promoted by S
 +
kIR=0.39*conv;%0.39*conv;    %Constant k of the hill equation for the promoter promoted by the complex luxIluxR
 +
kCI=0.055*conv;    %Constant k of the hill equation for the promoter promoted by CI
 +
%
 +
hS=1;    %Hill constant for the promoters promoted by S
 +
hIR=3.4;    %Hill constant for the promoter promoted by the complex IR
 +
hCI=2.3;    %Hill constant fot the promoter CI
 +
%
 +
eA=0.5/conv;    %Export factor of the chitinase
 +
jQ=0.1/conv^3;    %Import factor of the chitin monomers
 +
deltaA=0.2/conv; %Difusion factor of the chinitanse outside the cell
 +
%
 +
eI=0.5/conv;    %Export factor of LuxI
 +
jI=0.8/conv;    %Import factor of LuxI
 +
deltaI=0.2/conv;%Difusion of LuxI outside the cell
 +
%
 +
Stotal= 1.5*conv;  %Total concentration of the sensor in the cell
 +
%
 +
eAS=0.8/conv;
 +
%
 +
%
 +
%-------------------------------------------------------------------------%
 +
%                      Time and number of cells conditions
 +
%-------------------------------------------------------------------------%
 +
%
 +
numcel=3;
 +
numpas=2000000;
 +
%
 +
%-------------------------------------------------------------------------%
 +
%                              Initial conditions
 +
%-------------------------------------------------------------------------%
 +
%
 +
T=zeros(1,numpas);
 +
C=zeros(1,numpas);
 +
CS=zeros(1,numpas);
 +
P=zeros(1,numpas);
 +
Ai=zeros(1,numpas);
 +
Ao=zeros(1,numpas);
 +
Q=zeros(1,numpas);
 +
Ii=zeros(1,numpas);
 +
Io=zeros(1,numpas);
 +
IR=zeros(1,numpas);
 +
R=zeros(1,numpas);
 +
CI=zeros(1,numpas);
 +
HB=zeros(1,numpas);
 +
HA=zeros(1,numpas);
 +
AS=zeros(1,numpas);
 +
QQ=zeros(1,numpas);
 +
S=zeros(1,numpas);
 +
%
 +
%
 +
%
 +
C1=(-(gammaC + mCS*Stotal- alfC* mCS/gammaCS)+sqrt((gammaC + mCS*Stotal- alfC* mCS/gammaCS)^2+4*mCS*gammaC/gammaCS*alfC))/(2*mCS*gammaC/gammaCS);
 +
S1=Stotal/(1+mCS*C1/gammaCS);
 +
CS1=mCS*C1*S1/gammaCS; 
 +
P1=(alfP+ (betaP*(S1^hS))/(kS^hS+(S1^hS)))/gammaP;
 +
Ai1=(alfA+ (betaA*(S1^hS))/(kS^hS+(S1^hS)))/(eA+gammaA);
 +
Ao1= eA*Ai1/(deltaA);
 +
Q1= 0; %
 +
%
 +
xi=[1,1,1,(alfR/gammaR),(alfCI/gammaCI),1,1,(alfAS/gammaAS)];
 +
%
 +
yi=(fsolve(@CondIn2,xi));
 +
%
 +
C(1,1)=round(C1);%
 +
CS(1,1)=round(CS1);
 +
P(1,1)=round(P1);
 +
Ai(1,1)=round(Ai1);
 +
Ao(1,1)=round(Ao1);
 +
Q(1,1)=round(Q1);
 +
Ii(1,1)=round(abs(real(yi(1))));
 +
Io(1,1)=round(abs(real(yi(2))));
 +
IR(1,1)=round(abs(real(yi(3))));
 +
R(1,1)=round(abs(real(yi(4))));
 +
CI(1,1)=round(abs(real(yi(5))));
 +
HB(1,1)=round(abs(real(yi(6))));
 +
HA(1,1)=round(abs(real(yi(7))));
 +
AS(1,1)=round(abs(real(yi(8))));
 +
S(1,1)=Stotal - CS(1,1);%round( Stotal);
 +
%
 +
Tmaximo=1;
 +
%
 +
%-------------------------------------------------------------------------%
 +
%                              Operational Code
 +
%-------------------------------------------------------------------------%
 +
%
 +
for j=1:numcel
 +
% 
 +
    for i=1:numpas-1
 +
       
 +
        %-----------------------------------------------------------------%
 +
        %                          Chitin Pulses
 +
        %-----------------------------------------------------------------%
 +
       
 +
        QQ(1,i)=functionChi(T(1,i));
 +
       
 +
        %-----------------------------------------------------------------%
 +
        %                    EVENTS OF THE SYSTEM
 +
        %-----------------------------------------------------------------%
 +
       
 +
        a=alfC; %Creation of CBP
 +
       
 +
        b= gammaC*C(1,i); %Degradation of CBP
 +
       
 +
        c= mCS*C(1,i)*S(1,i); %Formation complex CS
 +
       
 +
        d= mCSQ*CS(1,i)*Q(1,i); %Reaction CS with Q
 +
       
 +
        e =gammaC*C(1,i);%Degradation of complex CS
 +
       
 +
        f=alfP  + (betaP*(S(1,i)^hS))/(kS^hS+(S(1,i)^hS));%Creation of chitoporin
 +
       
 +
        g=gammaP*P(1,i);  %Destruction of chitoporin
 +
       
 +
        h=alfA + (betaA*(S(1,i)^hS))/(kS^hS+(S(1,i)^hS)); %Creation chitinase inside the cell
 +
       
 +
        k=gammaA*Ai(1,i);%Destruction of chitinase inside the cell
 +
       
 +
        l=eA*Ai(1,i); %Exportation od chitinase
 +
       
 +
        m=deltaA*Ao(1,i); %Dilution outside the cell
 +
       
 +
        o=mAQQ*Ao(1,i)*QQ(1,i); %Reaction of chitinase with chitin
 +
       
 +
        p= 2*jQ*P(1,i)*(mAQQ*QQ(1,i)*Ao(1,i)); %Creation of chitin monomer inside the cell
 +
       
 +
        q= alfI+ (betaI*(S(1,i)^hS))/(kS^hS+(S(1,i)^hS));% Production of luxI
 +
       
 +
        r=jI*Io(1,i); %Import of lux I outside the cell
 +
       
 +
        s= gammaI*Ii(1,i); %Degradation of luxI
 +
       
 +
        t= eI*Ii(1,i); %Export of lux I
 +
       
 +
        u= mIR*Ii(1,i)*R(1,i); %Formation of the complex LuxI-LuxR
 +
       
 +
        v= deltaI*Io(1,i); %Dilution of LuxI outside the cell
 +
       
 +
        w= mI*IR(1,i);  %Destruction of the lux complex
 +
       
 +
        x= alfR+(betaI*(S(1,i)^hS))/(kS^hS+(S(1,i)^hS)); %Creation of LuxR
 +
       
 +
        y= gammaR*R(1,i); %Destruction of LuxR
 +
       
 +
        z= alfCI + (betaCI*(CI(1,i)^hCI))/(kCI^hCI+(CI(1,i)^hCI)) +(betaCI*(IR(1,i)^hIR))/(kIR^hIR+(IR(1,i)^hIR));%Creation of CI
 +
       
 +
        aa= gammaCI*CI(1,i);%Destruction of CI
 +
       
 +
        ab=alfHB +(betaHB*(CI(1,i)^hCI))/(kCI^hCI+(CI(1,i)^hCI))+(betaHB*(IR(1,i)^hIR))/(kIR^hIR+(IR(1,i)^hIR)); %Creation of HipB
 +
       
 +
        ac=gammaHB*HB(1,i); %Degradation of HipB
 +
       
 +
        if HA(1,i)>= 2 && HB(1,i)>= 2
 +
       
 +
        ad=mHAHB*HA(1,i)^2*HB(1,i)^2; %Destruction of HipB
 +
       
 +
        else
 +
           
 +
        ad=0;
 +
       
 +
        end
 +
       
 +
        af=alfHA+ (betaHA*(CI(1,i)^hCI))/(kCI^hCI+(CI(1,i)^hCI)); %Creation of HipA7
 +
       
 +
        ae=gammaHA*HA(1,i); %Destruction of HipA7
 +
       
 +
        ag=alfAS+(betaAS*(CI(1,i)^hCI))/(kCI^hCI+(CI(1,i)^hCI))+(betaHA*(CI(1,i)^hCI))/(kCI^hCI+(CI(1,i)^hCI));  %Creation of Salicylic acid
 +
       
 +
        ah= gammaAS*AS(1,i); %Degradation of Salicylic acid
 +
       
 +
        ai=eAS*AS(1,i);  %Export of Salicylic acid
 +
       
 +
        %-----------------------------------------------------------------%
 +
        %                        Gillespie Algorithm
 +
        %-----------------------------------------------------------------%
 +
       
 +
        med=a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad+ae+af+ag+ah+ai;
 +
       
 +
        T(1,i+1)=T(1,i)+1/(med) * log(1/rand());    %Time distribution
 +
       
 +
        n=rand();
 +
       
 +
       
 +
        if (n>0) && (n<a/med)      %A molecule of C is created
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i)+1;
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
           
 +
        elseif (n>a/med) && (n<(a+b)/med) && (C(1,i)>0)  %A molecule of C is consumed
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i)-1;
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b)/med) && (n<(a+b+c)/med) && C(1,i)>0 %The complex C-S is created
 +
           
 +
            S(1,i)=Stotal-1;
 +
            C(1,i+1)=C(1,i)-1;
 +
            CS(1,i+1)=CS(1,i)+1;
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c)/med) && (n<(a+b+c+d)/med) && (CS(1,i)>0) && QQ(1,i)>0 %The complex C-S attaches to chitin and frees the sensor
 +
           
 +
            S(1,i)=Stotal+1;
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i)-1;
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i)-1;
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d)/med) && (n<(a+b+c+d+e)/med) && CS(1,i)>0 && QQ(1,i)>0 %The complex CS is destroyed
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i)-1;
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e)/med) && (n<(a+b+c+d+e+f)/med)  %Chitoporin is created
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i)+1;
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f)/med) && (n<(a+b+c+d+e+f+g)/med) && P(1,i)>0 %The chitoporine is diluted
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i)-1;
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g)/med) && (n<(a+b+c+d+e+f+g+h)/med) % Thei chitinase is created
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i)+1;
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h)/med) && (n<(a+b+c+d+e+f+g+h+k)/med) && Ai(1,i)>0 %The chitinase inside the cell is destructed
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i)-1;
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k)/med) && (n<(a+b+c+d+e+f+g+h+k+l)/med) && (Ai(1,i)>0) % The chitinase is exported
 +
           
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i)-1;
 +
            Ao(1,i+1)=Ao(1,i)+1;
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k+l)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m)/med)&& Ao(1,i)>0 %the hitinase outside the cell is diluted
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i)-1;
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k+l+m)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o)/med) && (Ao(1,i)>0) && QQ(1,i)>0 %The chitinase clivajes the chitin
 +
           
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i)-1;
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p)/med) && QQ(1,i)>0 %Thi chitin enters to cell
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i)+1;
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q)/med)  %The LuxI concentration inside the cell increases
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i)+1;
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r)/med) && Io(1,i)>0 %Import of LuxI
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i)+1;
 +
            Io(1,i+1)=Io(1,i)-1;
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s)/med) && (Ii(1,i)>0)  %The LuxI concentration inside the cell decreases
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i)-1;
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t)/med) && Ii(1,i)>0 %Export of LuxI
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i)-1;
 +
            Io(1,i+1)=Io(1,i)+1;
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u)/med) && (Ii(1,i)>0) && R(1,i)>0 %The complex LuxI-LuR is created
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i)-1;
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i)+1;
 +
            R(1,i+1)=R(1,i)-1;
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v)/med) &&Io(1,i)>0 %Dilution of LuxI outsed the cell
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i)-1;
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w)/med) && (IR(1,i)>0)  %LuxI-LuxR is destryed
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i)-1;
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x)/med) % LuxR is created
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i)+1;
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y)/med) && (R(1,i)>0)  %R decresases
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i)-1;
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z)/med) %CI is created
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i)+1;
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa)/med) && (CI(1,i)>0)  %CI decreases
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i)-1;
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab)/med) %HB increases
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i)+1;
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac)/med) && (HB(1,i)>0)  %HB decreases
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i)-1;
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad)/med) && HA(1,i)>1 && HB(1,i)>1 %Inactivation of the toxin
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i)-2;
 +
            HA(1,i+1)=HA(1,i)-2;
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad+ae)/med) && (HA(1,i)>0)  %HA decreases
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i)-1;
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad+ae)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad+ae+af)/med)  %HA decreases
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i)+1;
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad+ae+af)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad+ae+af+ag)/med)  %AS Increases
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i)+1;
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad+ae+af+ag)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad+ae+af+ag+ah)/med) && AS(1,i)>0  %AS Decreases
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i)-1;
 +
           
 +
        elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad+ae+af+ag+ah)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad+ae+af+ag+ah+ai)/med) && AS(1,i)>0 %AS Decreases
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i)-1;
 +
           
 +
        else
 +
           
 +
            S(1,i+1)= S(1,i);
 +
            C(1,i+1)=C(1,i);
 +
            CS(1,i+1)=CS(1,i);
 +
            P(1,i+1)=P(1,i);
 +
            Ai(1,i+1)=Ai(1,i);
 +
            Ao(1,i+1)=Ao(1,i);
 +
            Q(1,i+1)=Q(1,i);
 +
            Ii(1,i+1)=Ii(1,i);
 +
            Io(1,i+1)=Io(1,i);
 +
            IR(1,i+1)=IR(1,i);
 +
            R(1,i+1)=R(1,i);
 +
            CI(1,i+1)=CI(1,i);
 +
            HB(1,i+1)=HB(1,i);
 +
            HA(1,i+1)=HA(1,i);
 +
            AS(1,i+1)=AS(1,i);
 +
           
 +
           
 +
           
 +
           
 +
           
 +
        end
 +
       
 +
        if T(1,i) >= Tmaximo
 +
           
 +
            break
 +
           
 +
        end
 +
       
 +
       
 +
    end
 +
   
 +
    %---------------------------------------------------------------------%
 +
    %                      Converting to regular range
 +
    %---------------------------------------------------------------------%
 +
   
 +
    tmax=Tmaximo;%max(T);
 +
    stepp=0.001;
 +
   
 +
    [TFixed, CFixed]=regintervalfixed(T,C,stepp,tmax);
 +
    [TFixed, CSFixed]=regintervalfixed(T,CS,stepp,tmax);
 +
    [TFixed, PFixed]=regintervalfixed(T,P,stepp,tmax);
 +
    [TFixed, AiFixed]=regintervalfixed(T,Ai,stepp,tmax);
 +
    [TFixed, AoFixed]=regintervalfixed(T,Ao,stepp,tmax);
 +
    [TFixed, QFixed]=regintervalfixed(T,Q,stepp,tmax);
 +
    [TFixed, IiFixed]=regintervalfixed(T,Ii,stepp,tmax);
 +
    [TFixed, IoFixed]=regintervalfixed(T,Io,stepp,tmax);
 +
    [TFixed, IRFixed]=regintervalfixed(T,IR,stepp,tmax);
 +
    [TFixed, RFixed]=regintervalfixed(T,R,stepp,tmax);
 +
    [TFixed, CIFixed]=regintervalfixed(T,CI,stepp,tmax);
 +
    [TFixed, HBFixed]=regintervalfixed(T,HB,stepp,tmax);
 +
    [TFixed, HAFixed]=regintervalfixed(T,HA,stepp,tmax);
 +
    [TFixed, ASFixed]=regintervalfixed(T,AS,stepp,tmax);
 +
    [TFixed, QQFixed]=regintervalfixed(T,QQ,stepp,tmax);
 +
    [TFixed, SFixed]=regintervalfixed(T,S,stepp,tmax);
 +
   
 +
    for k=1:length(TFixed)
 +
       
 +
       
 +
        TMatrix(j,k)=TFixed(1,k);
 +
        CMatrix(j,k)=CFixed(1,k);
 +
        CSMatrix(j,k)=CSFixed(1,k);
 +
        PMatrix(j,k)=PFixed(1,k);
 +
        AiMatrix(j,k)=AiFixed(1,k);
 +
        AoMatrix(j,k)=AoFixed(1,k);
 +
        QMatrix(j,k)=QFixed(1,k);
 +
        IiMatrix(j,k)=IiFixed(1,k);
 +
        IoMatrix(j,k)=IoFixed(1,k);
 +
        IRMatrix(j,k)=IRFixed(1,k);
 +
        RMatrix(j,k)=RFixed(1,k);
 +
        CIMatrix(j,k)=CIFixed(1,k);
 +
        HBMatrix(j,k)=HBFixed(1,k);
 +
        HAMatrix(j,k)=HAFixed(1,k);
 +
        ASMatrix(j,k)=ASFixed(1,k);
 +
        QQMatrix(j,k)=QQFixed(1,k);
 +
        SMatrix(j,k)=SFixed(1,k);
 +
       
 +
       
 +
    end
 +
   
 +
   
 +
end
 +
%
 +
%-------------------------------------------------------------------------%
 +
%                          Mean values for answers
 +
%-------------------------------------------------------------------------%
 +
%
 +
%
 +
%
 +
%
 +
%
 +
TPromedio=mean(TMatrix);
 +
CPromedio=mean(CMatrix);
 +
CSPromedio=mean(CSMatrix);
 +
PPromedio=mean(PMatrix);
 +
AiPromedio=mean(AiMatrix);
 +
AoPromedio=mean(AoMatrix);
 +
QPromedio=mean(QMatrix);
 +
IiPromedio=mean(IiMatrix);
 +
IoPromedio=mean(IoMatrix);
 +
IRPromedio=mean(IRMatrix);
 +
RPromedio=mean(RMatrix);
 +
CIPromedio=mean(CIMatrix);
 +
HBPromedio=mean(HBMatrix);
 +
HAPromedio=mean(HAMatrix);
 +
ASPromedio=mean(ASMatrix);
 +
QQPromedio=mean(QQMatrix);
 +
SPromedio=mean(SMatrix);
 +
%
 +
%-------------------------------------------------------------------------%
 +
%                      Standard deviation for answers
 +
%-------------------------------------------------------------------------%
 +
%
 +
CDesv=std(CMatrix);
 +
CSDesv=std(CSMatrix);
 +
PDesv=std(PMatrix);
 +
AiDesv=std(AiMatrix);
 +
AoDesv=std(AoMatrix);
 +
QDesv=std(QMatrix);
 +
IiDesv=std(IiMatrix);
 +
IoDesv=std(IoMatrix);
 +
IRDesv=std(IRMatrix);
 +
RDesv=std(RMatrix);
 +
CIDesv=std(CIMatrix);
 +
HBDesv=std(HBMatrix);
 +
HADesv=std(HAMatrix);
 +
ASDesv=std(ASMatrix);
 +
QQDesv=std(QQMatrix);
 +
SDesv=std(SMatrix);
 +
%
 +
%-------------------------------------------------------------------------%
 +
%                                Plotting
 +
%-------------------------------------------------------------------------%
 +
%
 +
plot(TPromedio,ASPromedio,TPromedio,CIPromedio,TPromedio,QQPromedio/50);
 +
xlabel('Time')
 +
%
 +
ylabel('Number of Molecules')
 +
%legend('CBP','CBP-S','Chitoporin','Chitinase inside','Chitinase outside','Chitin monomers','Quitin');
 +
%
 +
%
 +
%
 +
%
 +
%
 +
%
 +
%
 +
%

Revision as of 03:47, 27 October 2012

Template:Https://2012.igem.org/User:Tabima

Contents

Scripting

RALSTONIA DIFFERENTIAL EQUATION SOLUTION

These scripts create the differential equations, find the steady state concetrations and then solve them by a 4th order Runge Kutta:

DIFFERANTIAL EQUATIONS DECLARATIONS

%THIS CODE CREATE ALL THE DIFFERENTIAL EQUATIONS FOR RALSTONIA'S SYSTEM 

function y=ecuaDifR(t,v)                        
                                                        %---------Parameters------%
global alfS   %Basal concentration of the sensor phcS
global alfRA   %Basal concnetration of the comple pchA-pchR 
global alfR   %Basal concentration of LuxR
global alfI   %Basal concentration of CI 
global alfCI  %Basal concentration of CI 
global alfHA  %Basal concnetration of HipA7
global alfHB  %Basal concnetration of HipB 
global alfAS   %Basal concnetration of Salycilic acid 
global gammaS   %Degradation of the sensor pchS
global gammaRA   %Degradation of the complex pchR-pchA
global gammaR   %Degradation of LuxR
global gammaI   %Degradation of LuxI 
global gammaCI  %Degradation of CI 
global gammaHA  %Degradation of HipA7
global gammaHB  %Degradation of HipB 
global gammaAS  %Dergradation of Salycilic acid 
global mOHS     %Kinetic constant for the detection of 3-OH-PAME by the sensor pchS (phosphorylation)
global mSFR     %Kinetic constant for the phosphorylation of the complex pchR-pchA by the sensor
global mA       %Kinetic constant for the activation of the promoter by the pchA
global mIR      %Kinetic constant for the formation of the complex LuxILuxR
global mI       %Constant that represent the union of the complex LuxILuxR with the promoter
global mHAHB    %Kinetic constant for the inhibition of HipA7 
global betaI     %Max production of LuxI
global betaCI    %Max production of CI 
global betaHB    %Max peoduction of HipB
global betaHA    %Max production of HipA7 
global betaAS    %Max production of Salicylic acid
global kA       %Constant k of the hill ecuation for the promoter promoted by pchA
global kIR      %Constant k of the hill equiation for the promorer prmoted by the complex luxIluxR
global kCI      %Constant k of the hill equation for the promoter promoted by CI 
global hA      %Hill constant for the promoters promoted by pchA 
global hIR     %Hill constant for the promoter promoted by the complex IR
global hCI     %Hill constant fot the promoter CI 
global eI       %Export factor of LuxI
global jI       %Import factor of LuxI
global deltaI   %Difusion of LuxI outside the cell 
global eAS     %Export of Salicylic acid  
global numcel  %number of cells 

if (t<(10) || ((t)>20))
    
    OH=0;
    
else
    
    OH=15;
    
    
end                  

                    %------ Variables%------
                       
S=v(1); %Cocentration of the sensor pchS the cell
SF=v(2); %Concentration of phosphorylated sensor the cell
RA=v(3);  %Concentration of the comple pchR-pchA
A=v(4);  %Concentratio of the promoter avtivator pchA
Ii=v(5); %Concentration of LuxI inside the cell
Io=v(6); %Concentration of LuI outsied the cell
IR=v(7); %Concentration of the complex LuxI-LuxR 
R=v(8);%Concentration of the protein CI 
CI=v(9);%Concentration of HipA7
HB=v(10);%Concnetratio of HipB
HA=v(11);%Concentration of salicylic acid
AS=v(12); %Concentratio of quitin monomers
                                                 
                                %---Equations---%
             
                                
dS=alfS- gammaS*S - mOHS*OH*S; %Change of the sensor pchS
dSF = mOHS *OH*S - mSFR *SF*RA ;  %Change of phosphorylated sensor
dRA=alfRA - gammaRA*RA - mSFR*SF*RA;%Change of the comple pchR-pchA
dA= mSFR*SF*RA-mA*A; %Change of the activator pchA inside the cell
dIi= alfI+ (betaI*(A^hA))/(kA^hA+(A^hA)) -gammaI*Ii +jI*Io- eI*Ii- mIR*Ii*R; %Change of LuxI inside 
 the cell
dIo= numcel*(eI*Ii-jI*Io)-deltaI*Io; %Change of LuxI outside the cell
dIR= mIR*Ii*R - mI*IR;  %Change of the complex LuxI luxR
dR= alfR-gammaR*R -mIR*Ii*R +(betaI*(A^hA))/(kA^hA+(A^hA)); %Change of LuxR
dCI= alfCI -gammaCI*CI+ (betaCI*(CI^hCI))/(kCI^hCI+(CI^hCI)) +(betaCI*(IR^hIR))/(kIR^hIR+(IR^hIR));
%Change of CI
dHB=alfHB-gammaHB*HB+(betaHB*(CI^hCI))/(kCI^hCI+(CI^hCI))+(betaHB*(IR^hIR))/(kIR^hIR+(IR^hIR))
-mHAHB*HA^2*HB^2;  %Chanche of HipB 
dHA=alfHA-gammaHA*HA+ (betaHA*(CI^hCI))/(kCI^hCI+(CI^hCI))-mHAHB*HA^2*HB^2; %Change of HipA7
dAS=alfAS-gammaAS*AS +(betaAS*(CI^hCI))/(kCI^hCI+(CI^hCI))-eAS*AS+(betaAS*(IR^hIR))/(kIR^hIR+(IR^hIR));
%Change of Salicylic acid 
y1(1)=dS;
y1(2)=dSF;
y1(3)=dRA;
y1(4)=dA;
y1(5)=dIi;
y1(6)=dIo;
y1(7)=dIR;
y1(8)=dR;
y1(9)=dCI;
y1(10)=dHB;
y1(11)=dHA;
y1(12)=dAS; 
y=y1';         
                                               
end


RUNGE KUTTA 

%FILE THAT SOLVES THE DIFFERENTIAL EQUATION AND GRAPHS THEM 
alfS=0.9;      %Basal concentration of the sensor pchS
alfRA=0.9;      %Basal concnetration of the complez pchA-pchR
alfR=0.6;      %Basal concentration of LuxR
alfI=0.4;      %Basal concentration of LuxI
alfCI=0.5;     %Basal concentration of CI
alfHA=1;     %Basal concnetration of HipA7
alfHB=0.4;     %Basal concnetration of HipB
alfAS=0.4;    %Basal concnetration of Salycilic acid
gammaS=1;    %Degradation of Chitinase inside the cell
gammaRA=1;    %Degradation of chitoporin
gammaR=1;    %Degradation of LuxR
gammaI=1;    %Degradation of LuxI
gammaCI=1;   %Degradation of CI
gammaHA=1;   %Degradation of HipA7
gammaHB=4;   %Degradation of HipB
gammaAS=1;   %Dergradation of Salycilic acid 
mOHS=4;      %Kinetic constant for the formation of the phoshorilation of the sensor
mSFR=2.6;     %Kinetic constant of the reaction of the phosphorilarion of the comple R
mA=3.5;     %Kinetic constant for the activation by A 
mIR=3;      %Kinetic constant for the formation of the complex LuxILuxR
mI=3;       %Constant that represent the union of the complex LuxILuxR with the promoter
mHAHB=12;    %Kinetic constant for the inhibition of HipA7
betaI=10;    %Max production of LuxI
betaCI=9.96;   %Max production of CI
betaHB=9.95;   %Max production of HipB
betaHA=10;   %Max production of HipA7
betaAS=11.2;   %Max production of Salicylic acid
kA=0.1;      %Constant k of the hill ecuation for the promoter promoted by S
kIR=0.39;     %Constant k of the hill equiation for the promorer prmoted by the complex luxIluxR
kCI=0.055;     %Cosntant k of the hill equation for the promoter promoted by CI
hA=1.2;     %Hill constant for the promoters promoted by S
hIR=3.4;    %Hill constant for the promoter promoted by the complex IR
hCI=2.3;    %Hill constant fot the promoter CI
eI=0.5;    %Export factor of LuxI
jI=0.8;    %Import factor of LuxI
deltaI=0.2;%Difusion of LuxI outside the cell 
eAS=0.8;
numcel=1;   %number of cells
%----%

h=100; %Tiempo maximo

m=0.01; %Longitud de paso [s] 

t=0:m:h; %Vector tiempo

xi=[1,1,1,1,1,1,1,1,1,1,1,1];

y=fsolve(@CondInR,xi);

conInd=y;

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

OH=zeros(1,length(l));

x(1,:)=conInd;

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=ecuaDifR(l(k),xk); %Primera pendiente del metodo de RK4
   k2=ecuaDifR(l(k)+m/2,xk+(m/2*k1)'); %Segunda pendiente del metodo de RK4
   k3=ecuaDifR(l(k)+m/2,xk+(m/2*k2)'); %Tercera pendiente del metodo de RK4
   k4=ecuaDifR(l(k)+m,xk+(m*k3)'); %Cuarta pendiente del metodo de RK4
   
   xk1=xk+m/6*(k1+2*k2+2*k3+k4)'; %Calculo de nuevos valores para las
   %variables
   
   %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

for j=1:length(l)
    
    if (l(j)<(10) || l(j)>(20))
       
       OH(j)=0;
       
   else
       
       OH(j)=15;
       
       
   end
   
   
end

S=x(:,1);
SF=x(:,2);
RA=x(:,3);
A=x(:,4);
Ii=x(:,5);
Io=x(:,6);
IR=x(:,7);
R=x(:,8);
CI=x(:,9);
HB=x(:,10);
HA=x(:,11);
AS=x(:,12); 

figure(1) 
plot(l,S,l,SF)
legend('Sensor (pchS)',' Phosporilated pchS')
xlabel('Time')
ylabel('Concetratio (micromolar)')
title('Response of Sensor pchS')
figure(5)
plot(l,RA,l,A)
legend('Complex pchR-pchA','Activator pchA')
xlabel('Time')
ylabel('Concetration (micromolar)')
title('Activator response')
figure(2)
plot (l,R,l,Ii,l,Io,l,IR)
legend('LuxR','LuxI nside the cell','LuxI outside the cell','Complex(Lux-LuxR)')
xlabel('Time')
ylabel('Concetration (micromolar)')
title('LuxI-LuxR system response')
figure(3) 
plot (l,HA,l,HB) 
legend('Toxin HipA7','Antitoxin HipB')
xlabel('Time')
ylabel('Concetration (micromolar)')
title('Toxin-Antitoxin module')
figure(4) 
plot (l,OH,l,AS,l,CI) 
legend('3-OH-PAME','Salicylic acid','CI')
xlabel('Time')
ylabel('Concetration (micromolar)')
title('CI and Salicylic Acid response')

RUST DIFFERENTIAL EQUATION SOLUTION

Chitin impulse function: This function returns the chitin concentration depending on the time imput 

function answer = functionChi( t )
%
answer=0;
%
if t>10 && t<20
 %   
   answer=15;
  % 
end
%

end


Differential equation declaration: 

%THIS CODE CREATE ALL THE DIFFERENTIAL EQUATIONS FOR THE SYSTEM
%
function y=ecuaDif(t,v)                        
                                                        %---------Parameters------%
global alfA   %Basal concentration of Chitinase inside the cell (micromolar)
global alfP   %Basal concnetration of chitoporin 
global alfC   %Basal concentration of the CBP
global alfR   %Basal concentration of LuxR
global alfI   %Basal concentration of CI 
global alfCI  %Basal concentration of CI 
global alfHA  %Basal concnetration of HipA7
global alfHB  %Basal concnetration of HipB 
global alfAS   %Basal concnetration of Salycilic acid 
global alfCS
global gammaA   %Degradation of Chitinase inside the cell
global gammaP   %Degradation of chitoporin 
global gammaC   %Degradation concentration of the CBP
global gammaR   %Degradation of LuxR
global gammaI   %Degradation of LuxI 
global gammaCI  %Degradation of CI 
global gammaHA  %Degradation of HipA7
global gammaHB  %Degradation of HipB 
global gammaAS  %Dergradation of Salycilic acid
global gammaCS  %Degradation of the complex CS 
global mCS      %Kinetic constant for the formation of the complex CS 
global mCSQ     %Kinetic constant of the reaction of the complex CS with the chitin 
global mAQQ     %Kinetic constant for the reaction of the chitinase and th chitin
global mIR      %Kinetic constant for the formation of the complex LuxILuxR
global mI       %Constant that represent the union of the complex LuxILuxR with the promoter
global mHAHB    %Kinetic constant for the inhibition of HipA7 
global betaP     %Max production of the chitoporin 
global betaA     %Max production of chitinase
global betaI     %Max production of LuxI
global betaCI    %Max production of CI 
global betaHB    %Max peoduction of HipB
global betaHA    %Max production of HipA7 
global betaAS    %Max production of Salicylic acid 
global kS       %Constant k of the hill ecuation for the promoter promoted by S
global kIR      %Constant k of the hill equiation for the promorer prmoted by the complex luxIluxR
global kCI      %Constant k of the hill equation for the promoter promoted by CI 
global hS      %Hill constant for the promoters promoted by S 
global hIR     %Hill constant for the promoter promoted by the complex IR
global hCI     %Hill constant fot the promoter CI  
global eA      %Export factor of the chitinase 
global jQ      %Import factor of the chitin monomers 
global deltaA  %Difusion factor of the chinitanse outside the cell
global eI       %Export factor of LuxI
global jI       %Import factor of LuxI
global deltaI   %Difusion of LuxI outside the cell  
global Stotal  %Total concentration of the sensor in the cell 
global eAS     %Export of Salicylic acid 
global numcel  %number of cells 
% 
% 
QQ=functionChi(t);
%
                           %------ Variables%------
%
%                       
C=v(1); %Cocentration of chitinase outside the cell
CS=v(2); %Concentration of chitinase inside the cell
P=v(3);  %Concentration of chitiporin
Ai=v(4);  %Concentratio of chitin binding protein (CBP)
Ao=v(5); %Concentration the complex CBP-s
Q=v(6);  %Concentration of LuxR
Ii=v(7); %Concentration of LuxI inside the cell
Io=v(8); %Concentration of LuI outsied the cell
IR=v(9); %Concentration of the complex LuxI-LuxR 
R=v(10);%Concentration of the protein CI 
CI=v(11);%Concentration of HipA7
HB=v(12);%Concnetratio of HipB
HA=v(13);%Concentration of salicylic acid
AS=v(14); %Concentratio of quitin monomers
%
%                                                
%                                %---Equations---%
                                
  $                            
  % 
S=Stotal-CS;
  %                             
  %                             
  %
dC=alfC- gammaC*C - mCS*C*S; %Change of CBP 
dCS=alfCS+ mCS*C*S- mCSQ*CS*Q-gammaCS*CS;  %Change of the complex CS
dP=alfP - gammaP*P + (betaP*(S^hS))/(kS^hS+(S^hS));%Change of chitoporin
dAi=(alfA- gammaA*Ai+ (betaA*(S^hS))/(kS^hS+(S^hS)))- eA*Ai; %Change of chitinase inside the cell 
dAo= eA*Ai-deltaA*Ao- mAQQ*Ao*QQ; %Change of chitinase outside the cell
dQ= 2*jQ*P*(mAQQ*QQ*Ao)-mCSQ*CS*Q; %Change of chitin monomer inside the cell
dIi= alfI+ (betaI*(S^hS))/(kS^hS+(S^hS)) -gammaI*Ii +jI*Io- eI*Ii- mIR*Ii*R; %Change of LuxI inside 
the cell
dIo= numcel*(eI*Ii-jI*Io)-deltaI*Io; %Change of LuxI outside the cell
dIR= mIR*Ii*R - mI*IR;  %Change of the complex LuxI luxR
dR= alfR-gammaR*R -mIR*Ii*R +(betaI*(S^hS))/(kS^hS+(S^hS)); %Change of LuxR
dCI= alfCI -gammaCI*CI+ (betaCI*(CI^hCI))/(kCI^hCI+(CI^hCI)) +(betaCI*(IR^hIR))/(kIR^hIR+(IR^hIR));
%Change of CI
dHB=alfHB-gammaHB*HB+(betaHB*(CI^hCI))/(kCI^hCI+(CI^hCI))+(betaHB*(IR^hIR))/(kIR^hIR+(IR^hIR))
-mHAHB*HA^2*HB^2;  %Chanche of HipB 
dHA=alfHA-gammaHA*HA-mHAHB*HA^2*HB^2+(betaHA*(CI^hCI))/(kCI^hCI+(CI^hCI)); %Change of HipA7
dAS=alfAS-gammaAS*AS +(betaCI*(CI^hCI))/(kCI^hCI+(CI^hCI))+(betaAS*(IR^hIR))/(kIR^hIR+(IR^hIR))-eAS*AS;
%Change of Salicylic acid
%
y1(1)=dC;
y1(2)=dCS;
y1(3)=dP;
y1(4)=dAi;
y1(5)=dAo;
y1(6)=dQ;
y1(7)=dIi;
y1(8)=dIo;
y1(9)=dIR;
y1(10)=dR;
y1(11)=dCI;
y1(12)=dHB;
y1(13)=dHA;
y1(14)=dAS;
%
y=y1';            
%                        
%                          
end


Diferential equation solution 

tic;
%
%File that solves the differential equations and graphs them
%
alfA = 0.9;      %Basal concentration of Chitinase inside the cell (micromolar)
alfP = 0.9;      %Basal concnetration of chitoporin
alfC = 0.9;      %Basal concentration of the CBP
alfR = 0.6;      %Basal concentration of LuxR
alfI = 0.4;      %Basal concentration of LuxI
alfCI = 0.5;     %Basal concentration of CI
alfHA = 1.4;     %Basal concnetration of HipA7
alfHB = 0.4;     %Basal concnetration of HipB
alfAS = 0.4;    %Basal concnetration of Salycilic acid
alfCS=1.4;
gammaA=1;    %Degradation of Chitinase inside the cell
gammaP=1;    %Degradation of chitoporin
gammaC=1;    %Degradation concentration of the CBP
gammaR=1;    %Degradation of LuxR
gammaI=1;    %Degradation of LuxI
gammaCI=1;   %Degradation of CI
gammaHA=1;   %Degradation of HipA7
gammaHB=4;   %Degradation of HipB
gammaAS=0.8;   %Dergradation of Salycilic acid
gammaCS=1;   %Degradation of the complex Cs 
mCS=13;      %Kinetic constant for the formation of the complex CS
mCSQ=12;     %Kinetic constant of the reaction of the complex CS with the chitin
mAQQ=0.2;     %Kinetic constant for the reaction of the chitinase and th chitin
mIR=3;      %Kinetic constant for the formation of the complex LuxILuxR
mI=3;       %Constant that represent the union of the complex LuxILuxR with the promoter
mHAHB=12;    %Kinetic constant for the inhibition of HipA7 
betaP=12;    %Max production of the chitoporin
betaA=10;    %Max production of chitinase
betaI=10;    %Max production of LuxI
betaCI=9.96;   %Max production of CI
betaHB=9.5;   %Max production of HipB
betaHA=11;   %Max production of HipA7
betaAS=11.2;   %Max production of Salicylic acid 
kS=0.08;      %Constant k of the hill ecuation for the promoter promoted by S
kIR=0.39;     %Constant k of the hill equiation for the promorer prmoted by the complex luxIluxR
kCI=0.055;     %Cosntant k of the hill equation for the promoter promoted by CI 
hS=1;     %Hill constant for the promoters promoted by S
hIR=3.4;    %Hill constant for the promoter promoted by the complex IR
hCI=2.3;    %Hill constant fot the promoter CI
eA=0.5;     %Export factor of the chitinase
jQ=0.1;     %Import factor of the chitin monomers
deltaA=0.2; %Difusion factor of the chinitanse outside the cell 
eI=0.5;    %Export factor of LuxI
jI=0.8;    %Import factor of LuxI
deltaI=0.2;%Difusion of LuxI outside the cell
Stotal= 1.5;  %Total concentration of the sensor in the cell
eAS=0.8;
numcel=1;   %number of cells
%
h=100; %Tiempo maximo
%
m=0.01; %Longitud de paso [s]
%
t=0:m:h; %Vector tiempo
%
xi=[1,1,1,1,1,1,1,1,1,1,1,1,1,1];
%
y=fsolve(@CondIn,xi);
%
%
conInd=y;
%
%
%
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 la longitud
%
QQ=zeros(1,length(l));
%
x1=conInd;
%
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=ecuaDif(l(k),xk); %Primera pendiente del metodo de RK4
   k2=ecuaDif(l(k)+m/2,xk+(m/2*k1)'); %Segunda pendiente del metodo de RK4
   k3=ecuaDif(l(k)+m/2,xk+(m/2*k2)'); %Tercera pendiente del metodo de RK4
   k4=ecuaDif(l(k)+m,xk+(m*k3)'); %Cuarta pendiente del metodo de RK4
   
   xk1=xk+m/6*(k1+2*k2+2*k3+k4)'; %Calculo de nuevos valores para las
   %variables
   
   %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
% 
%
%
for j=1:length(l)
 %   
   QQ(1,j)=functionChi(l(j));
   % 
end
%
%
C=x(:,1);
CS=x(:,2);
P=x(:,3);
Ai=x(:,4);
Ao=x(:,5);
Q=x(:,6);
Ii=x(:,7);
Io=x(:,8);
IR=x(:,9);
R=x(:,10);
CI=x(:,11);
HB=x(:,12);
HA=x(:,13);
AS=x(:,14); 
%
%
figure(1) 
plot(l,C,l,CS,l,P,l,Ai,l,Ao,l,Q,l,QQ)
legend('CBP','Complex Sensor(CBP)','Chitoporin','Chitinase inside the cell','Chitinase outside the
cell','Chitin   monomers','Chitin')
xlabel('Time')
ylabel('Concentration')
title('Detection system substances')
%
figure(2)
plot (l,R,l,Ii,l,Io,l,IR)
legend('LuxR','LuxI Inside the cell','LuxI outside the cell','Complex LuxI-LuxR')
xlabel('Time')
ylabel('Concentration')
title('LuxI-LuxR system substances') 
%
figure(3) 
plot (l,HA,l,HB) 
legend('HipA7','HipB')
xlabel('Time')
ylabel('Concentration')
title('Toxin-Antitoxin module substances')
%
figure(4) 
plot (l,QQ,l,AS,l ,CI) 
legend('QQ','Salicylic acid','CI')
xlabel('Time')
ylabel('Concentration')
title('CI and Salicylic acid response')
%

OPTIMIZATION

The optimization was done using the software GAMS. Here we present the code: 

sets
t /1*201/
;
*
variable
z       Objective function
;
*
positive variables
alfc    constituive production of CBP
alfp    constitutive production chitoporin
alfa    constitutive production chitinase
alfi    basal production luxI
alfr    basal production luxR
alfci   basal production CI
alfha  basal production HipA7
alfhb   basal production HipB
alfas   basal production Salicylic acid
mcs     Kinetic constant
mcsq    Kinetic constant
maqq    Kinetic constant
mir     Kinetic constant
mi      Kinetic constant
mha     Kinetic constant
betap   Maximal production of chitoporin
betaa   Maximal production of chitinase
betai   Maximal production of LuxI
betahb  Maximal production of HipB
betaha  Maximal production of HipA7
betaas  Maximal production of Salicylic acid
ks      Hill's constant for the activation by S
kir     Hill's constant for the activation by IR
hs      Hill's constant for the activation by S
hir     Hill's constant for the activation by IR
ea      Export factor of chitinase
jq      Import factor of chitine
deltaa  Diffusion factor of chitinase
ei      Export factor of LuxI
ji      Import factor of LuxI
deltai  Difussion factor of LuxI
stotal  Total concetration of the sensor in the cell
eas     Export factor os Salicylic acid 
*
c(t)    Concentration of CBP
cs(t)   Concentration of complex CBP-S
p(t)    Concentration of chitoporin
ai(t)   Concentration of chitinase
ao(t)   Concentration of chitinase outside the cell
q(t)    Concentration of chitine monomer
ii(t)   Concentration of LuxI inside the cell
io(t)   Concentration of LuxI outside the cell
ir(t)   Concentration of complex LuxI-LuxR
r(t)    Concentration of LuxR
ci(t)   Concentration of CI
hb(t)   Concentration of HB
ha(t)   Concentration of HA
as(t)   Concentration of Salicylic Acid
qq(t)   Concentration of Chitin
s(t) 
;
*
alfc.lo=0.7;
alfc.up=1.2;
alfc.l=1;
*
alfp.lo=0.7;
alfp.up=1.2;
alfp.l=1;
*
alfa.lo=0.7;
alfa.up=1.2;
alfa.l=1;
*
alfi.lo=0.01;
alfi.up=0.6;
alfi.l=0.4;
*
alfr.lo=0.01;
alfr.up=0.6;
alfr.l=0.4;
*
alfci.lo=0.01;
alfci.up=0.6;
alfci.l=0.4;
*
alfha.lo=0.01;
alfha.up=0.6;
alfha.l=0.4;
*
alfhb.lo=0.01;
alfhb.up=0.6;
alfhb.l=0.4;
*
alfas.lo=0.01;
alfas.up=0.6;
alfas.l=0.4; 
*
mcs.lo=50;
mcs.up=1000;
mcs.l=500;
*
mcsq.up=1000;
mcsq.l=500;
*
maqq.up=1000;
maqq.l=500;
*
mir.up=1000;
mir.l=500;
*
mi.up=1000;
mi.l=20; 
*
mha.up=1000;
mha.l=500; 
*
betap.lo=5;
betap.up=23;
betap.l=14;
*
betaa.lo=5;
betaa.up=23;
betaa.l=14;
*
betai.lo=5;
betai.up=23;
betai.l=14;
*
betahb.lo=5;
betahb.up=23;
betahb.l=14;
*
betaha.lo=5;
betaha.up=23;
betaha.l=14; 
*
betaas.lo=5;
betaas.up=23;
betaas.l=14;
* 
ks.lo=0.01;
ks.up=0.9;
ks.l=0.05;
*
kir.lo=0.01;
kir.up=0.9;
kir.l=0.05;
*
hs.up=3;
hs.l=1;
*
hir.up=5;
hir.l=3; 
*
ea.lo=0.01;
ea.up=1;
ea.l=0.05;
*
jq.lo=0.01;
jq.up=1;
jq.l=0.05;
*
deltaa.lo=0.01;
deltaa.up=1;
deltaa.l=0.05;
*
ei.lo=0.01;
ei.up=1;
ei.l=0.05;
*
ji.lo=0.01;
ji.up=1;
ji.l=0.05;
*
deltai.lo=0.01;
deltai.up=1; 
deltai.l=0.05;
*
eas.lo=0.01;
eas.up=1;
eas.l=0.05;
*
stotal.lo=0.4;
stotal.up=2.5;
stotal.l=1;
*
c.l(t)=0.1616;
cs.l(t)=0.738;
p.l(t)=4.7233;
ai.l(t)=3.1489;
ao.l(t)=7.8722;
q.l(t)=0;
ii.l(t)=0.9662;
io.l(t)=0.4831;
ir.l(t)=3.6605;
r.l(t)=1.2628;
ci.l(t)=20.2;
hb.l(t)=0.5722;
ha.l(t)=2.5117;
as.l(t)=0.16933;
s.l(t)=1;
*
*
scalar
*
gammaa /1/
gammac /1/
gammap /1/
gammar /1/
gammai /1/
gammaci /1/
gammaha /1/
gammahb /4/
gammaas /1/
gammacs /1/ 
*
betaci /9.96/
kci /0.055/
hci /2.3/
*
numcel /1/
*
dt /0.5/
*
*
qq1 /0/
qq2 /20/
qq3 /60/
qq4 /5/ 
;
*
equations
*
cbss  Steady state equation of CBP
css
pss
aiss
aoss
qss
iiss
ioss
irss
rss
ciss
hbss
hass
asss
dc(t) Diferential equations
dcs(t)
dp(t)
dai(t)
dii(t)
dio(t)
dir(t)
dr(t)
dci(t)
dhb(t)
dha(t)
das(t)
dao1(t)
dq1(t)
dao2(t)
dq2(t)
dao3(t)
dq3(t)
sf(t)    function of the change of S in the cell
*
res1     restricion 1
res2     restriction 2
res3     restriction 3
*
*
fobj     Objective function
;
*
*Steady state equations
*
cbss.. alfc-gammac*c('1')-mcs*c('1')*s('1')=e=0;
css..  mcs*c('1')*s('1')-mcsq*cs('1')*q('1')-gammacs*cs('1')=e=0;
pss..  alfp-gammap*p('1')+ betap*(s('1')**hs)/(ks**hs+s('1')**hs)=e=0;
aiss.. alfa- gammaa*ai('1')+ betaa* (s('1')**hs)/((ks**hs)+(s('1')**hs))-ea*ai('1')=e=0;
aoss.. ea*ai('1')-deltaa*ao('1')=e=0;
qss.. -mcsq*cs('1')*q('1')=e=0;
iiss.. alfi + betai* (s('1')**hs)/((ks**hs)+ (s('1')**hs))-gammai*ii('1')+ji*io('1')-ei*ii('1')
-mir*ii('1')*r('1')=e=0;
ioss.. numcel*(ei*ii('1')-ji*io('1'))-deltai*io('1')=e=0;
irss.. mir*ii('1')*r('1')-mi*ir('1')=e=0;
rss.. alfr-gammar*r('1')- mir*ii('1')*r('1')+ betai* (s('1')**hs)/((ks**hs)+(s('1')**hs))=e=0;
ciss.. alfci-gammaci*ci('1')+ betaci* (ci('1')**hci)/((kci**hci)+(ci('1')**hci))+
betaci* (ir('1')**hir)/((kir**hir)+(ir('1')**hir))=e=0;
hbss.. alfhb-gammahb*hb('1')+ betahb* (ci('1')**hci)/((kci**hci)+(ci('1')**hci))+ 
betahb* (ir('1')**hir)/((kir**hir)+ (ir('1')**hir))-mha*power(ha('1'),2)*power(hb('1'),2)=e=0;
hass.. alfha-gammaha*ha('1')+ betaha* (ci('1')**hci)/((kci**hci)+(ci('1')**hci))
-mha*power(ha('1'),2)*power(hb('1'),2)=e=0;
asss.. alfas-gammaas*as('1')+ betaas* (ci('1')**hci)/((kci**hci)+(ci('1')**hci))- eas*as('1')=e=0;
* 
*Differential equations
*
dc(t).. alfc-gammac*c(t)-mcs*c(t)*s(t)=e=(c(t+1)-c(t))/dt;
dcs(t)..  mcs*c(t)*s(t)-mcsq*cs(t)*q(t)-gammacs*cs('1')=e=(cs(t+1)-cs(t))/dt ;
dp(t)..  alfp-gammap*p(t)+ betap* (s(t)**hs)/((ks**hs)+(s(t)**hs))=e=(p(t+1)-p(t))/dt;
dai(t).. alfa- gammaa*ai(t)+ betaa* (s(t)**hs)/((ks**hs)+(s(t)**hs))-ea*ai(t)=e=(ai(t+1)-ai(t))/dt ;
dii(t).. alfi + betai*(s(t)**hs)/((ks**hs)+(s(t)**hs))-gammai*ii(t)+ji*io(t)-ei*ii(t)
-mir*ii(t)*r(t)=e=(ii(t+1)-ii(t))/dt;
dio(t).. numcel*(ei*ii(t)-ji*io(t))-deltai*io(t)=e=(io(t+1)-io(t))/dt;
dir(t).. mir*ii(t)*r(t)-mi*ir(t)=e=(ir(t+1)-ir(t))/dt;
dr(t).. alfr-gammar*r(t)- mir*ii(t)*r(t)+ betai* (s(t)**hs)/((ks**hs)+(s(t)**hs))=e=(r(t+1)-r(t))/dt ;
dci(t).. alfci-gammaci*ci(t)+ betaci* (ci(t)**hci)/((kci**hci)+(ci(t)**hci))
+ betaci* (ir(t)**hir)/((kir**hir)+ (ir(t)**hir))=e=(ci(t+1)-ci(t))/dt ;
dhb(t).. alfhb-gammahb*hb(t)+ betahb* (ci(t)**hci)/((kci**hci)+(ci(t)**hci))
+ betahb* (ir(t)**hir)/((kir**hir)+(ir(t)**hir))-mha*power(ha(t),2)*power(hb(t),2)=e=(hb(t+1)-hb(t))/dt ;
dha(t).. alfha-gammaha*ha(t)+ betaha* (ci(t)**hci)/((kci**hci)+(ci(t)**hci))
-mha*power(ha(t),2)*power(hb(t),2)=e=(ha(t+1)-ha(t))/dt;
das(t).. alfas-gammaas*as(t)+ betaas* (ci(t)**hci)/((kci**hci)+(ci(t)**hci))
- eas*as(t)=e=(as(t+1)-as(t))/dt ;
*
*Differential equations depending on chitin impulse
*
dao1(t) $(ord(t)<51).. ea*ai(t)-deltaa*ao(t)-maqq*ao(t)*qq1=e=(ao(t+1)-ao(t))/dt;
dao2(t)$(ord(t)>50 and ord(t)<151).. ea* ai(t) -deltaa*ao(t)-maqq*ao(t)*qq2=e=(ao(t+1)-ao(t))/dt;
dao3(t) $(ord(t)>151 and ord(t)<201).. ea*ai(t)-deltaa*ao(t)-maqq*ao(t)*qq1=e=(ao(t+1)-ao(t))/dt;
dq1(t)$(ord(t)<51).. 2*jq*p(t)*maqq*qq1*ao(t) -mcsq*cs('1')*q('1')=e=(q(t+1)-q(t))/dt;
dq2(t)$(ord(t)>51 and ord(t)<151).. 2*jq*p(t)*maqq*qq2*ao(t) -mcsq*cs('1')*q('1')=e=(q(t+1)-q(t))/dt;
dq3(t)$(ord(t)>151 and ord(t)<201).. 2*jq*p(t)*maqq*qq1*ao(t) -mcsq*cs('1')*q('1')=e=(q(t+1)-q(t))/dt;
*S equation
sf(t).. s(t)=e=stotal-cs(t);
*Some restrictions
res1.. ha('1')=g= hb('1');
res2.. hb('131')=g=ha('131');
res3.. ha('181')=g=hb('181');
*objetive function
fobj.. z=e=power((as('1')-as('41')),2)+power((1.6*as('1')-as('111')),2)
+  power((1.8*as('1')-as('121')),2)+power((1.1*as('1')-as('161')),2)+power((as('1')-as('181')),2) ;
*
option nlp=ipopt;
model igem /all/;
solve igem using nlp minimizing z;


STOCHASTIC MODEL 

% MATLAB Code for modeling a gene net using Gillespie algorithm
% Code by Roberto Moran, Daniela Olivera, Cesar Quintana & Andrés Simbaqueba
% iGEM Colombia Team 2012, Universidad de Los Andes, Bogotá (Colombia)
%
%-------------------------------------------------------------------------%
%                                   Parameters
%-------------------------------------------------------------------------%
clear all;
%
conv=602.2; %Conversion factor from micromolar to molecules/cell
%
%
alfA = 0.9*conv;      %Basal production of Chitinase inside the cell (micromolar)
alfP = 0.9*conv;      %Basal production of chitoporin
alfC = 0.9*conv;      %Basal production of the CBP
alfR = 0.06*conv;      %Basal production of LuxR
alfI = 0.04*conv;      %Basal production of LuxI
alfCI = 0.01*conv;     %Basal production of CI
alfHA = 0.7*conv;     %Basal production of HipA7
alfHB = 0.2*conv;     %Basal production of HipB
alfAS = 0.4*conv;    %Basal production of Salycilic acid
%
gammaA=1;    %Degradation of Chitinase inside the cell
gammaP=1;    %Degradation of chitoporin
gammaC=1;    %Degradation concentration of the CBP
gammaR=1;    %Degradation of LuxR
gammaI=1;    %Degradation of LuxI
gammaCI=4;   %Degradation of CI
gammaHA=1;   %Degradation of HipA7
gammaHB=4;   %Degradation of HipB
gammaAS=0.8;   %Degradation of Salycilic acid
gammaCS=1;   %Degradation of the complex Cs
%
mCS=13/conv^2;      %Kinetic constant for the formation of the complex CS
mCSQ=12/conv^2;     %Kinetic constant of the reaction of the complex CS with the chitin
mAQQ=0.2/conv^2;     %Kinetic constant for the reaction of the chitinase and th chitin
mIR=3/conv^2;      %Kinetic constant for the formation of the complex LuxILuxR
mI=3/conv;       %Constant that represent the union of the complex LuxILuxR with the promoter
mHAHB=12/conv^4;    %Kinetic constant for the inhibition of HipA7
%
betaP=12*conv;    %Max production of the chitoporin
betaA=10*conv;    %Max production of chitinase
betaI=1*conv;    %Max production of LuxI
betaCI=1*conv;%9.96*conv;   %Max production of CI
betaHB=9.5*conv;   %Max production of HipB
betaHA=11*conv;   %Max production of HipA7
betaAS=1.12*conv;%11.2*conv;   %Max production of Salicylic acid
%
kS=0.08*conv;      %Constant k of the hill ecuation for the promoter promoted by S
kIR=0.39*conv;%0.39*conv;     %Constant k of the hill equation for the promoter promoted by the complex luxIluxR
kCI=0.055*conv;     %Constant k of the hill equation for the promoter promoted by CI
%
hS=1;     %Hill constant for the promoters promoted by S
hIR=3.4;    %Hill constant for the promoter promoted by the complex IR
hCI=2.3;    %Hill constant fot the promoter CI
%
eA=0.5/conv;     %Export factor of the chitinase
jQ=0.1/conv^3;     %Import factor of the chitin monomers
deltaA=0.2/conv; %Difusion factor of the chinitanse outside the cell
%
eI=0.5/conv;    %Export factor of LuxI
jI=0.8/conv;    %Import factor of LuxI
deltaI=0.2/conv;%Difusion of LuxI outside the cell
%
Stotal= 1.5*conv;  %Total concentration of the sensor in the cell
%
eAS=0.8/conv;
%
%
%-------------------------------------------------------------------------%
%                       Time and number of cells conditions
%-------------------------------------------------------------------------%
%
numcel=3;
numpas=2000000;
%
%-------------------------------------------------------------------------%
%                               Initial conditions
%-------------------------------------------------------------------------%
%
T=zeros(1,numpas);
C=zeros(1,numpas);
CS=zeros(1,numpas);
P=zeros(1,numpas);
Ai=zeros(1,numpas);
Ao=zeros(1,numpas);
Q=zeros(1,numpas);
Ii=zeros(1,numpas);
Io=zeros(1,numpas);
IR=zeros(1,numpas);
R=zeros(1,numpas);
CI=zeros(1,numpas);
HB=zeros(1,numpas);
HA=zeros(1,numpas);
AS=zeros(1,numpas);
QQ=zeros(1,numpas);
S=zeros(1,numpas);
%
%
%
C1=(-(gammaC + mCS*Stotal- alfC* mCS/gammaCS)+sqrt((gammaC + mCS*Stotal- alfC* mCS/gammaCS)^2+4*mCS*gammaC/gammaCS*alfC))/(2*mCS*gammaC/gammaCS);
S1=Stotal/(1+mCS*C1/gammaCS);
CS1=mCS*C1*S1/gammaCS;  
P1=(alfP+ (betaP*(S1^hS))/(kS^hS+(S1^hS)))/gammaP;
Ai1=(alfA+ (betaA*(S1^hS))/(kS^hS+(S1^hS)))/(eA+gammaA); 
Ao1= eA*Ai1/(deltaA);
Q1= 0; % 
%
xi=[1,1,1,(alfR/gammaR),(alfCI/gammaCI),1,1,(alfAS/gammaAS)];
%
yi=(fsolve(@CondIn2,xi));
%
C(1,1)=round(C1);%
CS(1,1)=round(CS1);
P(1,1)=round(P1);
Ai(1,1)=round(Ai1);
Ao(1,1)=round(Ao1);
Q(1,1)=round(Q1);
Ii(1,1)=round(abs(real(yi(1))));
Io(1,1)=round(abs(real(yi(2))));
IR(1,1)=round(abs(real(yi(3))));
R(1,1)=round(abs(real(yi(4))));
CI(1,1)=round(abs(real(yi(5))));
HB(1,1)=round(abs(real(yi(6))));
HA(1,1)=round(abs(real(yi(7))));
AS(1,1)=round(abs(real(yi(8))));
S(1,1)=Stotal - CS(1,1);%round( Stotal);
%
Tmaximo=1;
%
%-------------------------------------------------------------------------%
%                              Operational Code
%-------------------------------------------------------------------------%
%
for j=1:numcel
%   
   for i=1:numpas-1
       
       %-----------------------------------------------------------------%
       %                          Chitin Pulses
       %-----------------------------------------------------------------%
       
       QQ(1,i)=functionChi(T(1,i));
       
       %-----------------------------------------------------------------%
       %                     EVENTS OF THE SYSTEM
       %-----------------------------------------------------------------%
       
       a=alfC; %Creation of CBP
       
       b= gammaC*C(1,i); %Degradation of CBP
       
       c= mCS*C(1,i)*S(1,i); %Formation complex CS
       
       d= mCSQ*CS(1,i)*Q(1,i); %Reaction CS with Q
       
       e =gammaC*C(1,i);%Degradation of complex CS
       
       f=alfP  + (betaP*(S(1,i)^hS))/(kS^hS+(S(1,i)^hS));%Creation of chitoporin
       
       g=gammaP*P(1,i);  %Destruction of chitoporin
       
       h=alfA + (betaA*(S(1,i)^hS))/(kS^hS+(S(1,i)^hS)); %Creation chitinase inside the cell
       
       k=gammaA*Ai(1,i);%Destruction of chitinase inside the cell
       
       l=eA*Ai(1,i); %Exportation od chitinase
       
       m=deltaA*Ao(1,i); %Dilution outside the cell
       
       o=mAQQ*Ao(1,i)*QQ(1,i); %Reaction of chitinase with chitin
       
       p= 2*jQ*P(1,i)*(mAQQ*QQ(1,i)*Ao(1,i)); %Creation of chitin monomer inside the cell
       
       q= alfI+ (betaI*(S(1,i)^hS))/(kS^hS+(S(1,i)^hS));% Production of luxI
       
       r=jI*Io(1,i); %Import of lux I outside the cell
       
       s= gammaI*Ii(1,i); %Degradation of luxI
       
       t= eI*Ii(1,i); %Export of lux I
       
       u= mIR*Ii(1,i)*R(1,i); %Formation of the complex LuxI-LuxR
       
       v= deltaI*Io(1,i); %Dilution of LuxI outside the cell
       
       w= mI*IR(1,i);  %Destruction of the lux complex
       
       x= alfR+(betaI*(S(1,i)^hS))/(kS^hS+(S(1,i)^hS)); %Creation of LuxR
       
       y= gammaR*R(1,i); %Destruction of LuxR
       
       z= alfCI + (betaCI*(CI(1,i)^hCI))/(kCI^hCI+(CI(1,i)^hCI)) +(betaCI*(IR(1,i)^hIR))/(kIR^hIR+(IR(1,i)^hIR));%Creation of CI
       
       aa= gammaCI*CI(1,i);%Destruction of CI
       
       ab=alfHB +(betaHB*(CI(1,i)^hCI))/(kCI^hCI+(CI(1,i)^hCI))+(betaHB*(IR(1,i)^hIR))/(kIR^hIR+(IR(1,i)^hIR)); %Creation of HipB
       
       ac=gammaHB*HB(1,i); %Degradation of HipB
       
       if HA(1,i)>= 2 && HB(1,i)>= 2
       
       ad=mHAHB*HA(1,i)^2*HB(1,i)^2; %Destruction of HipB
       
       else
           
       ad=0;
       
       end
       
       af=alfHA+ (betaHA*(CI(1,i)^hCI))/(kCI^hCI+(CI(1,i)^hCI)); %Creation of HipA7
       
       ae=gammaHA*HA(1,i); %Destruction of HipA7
       
       ag=alfAS+(betaAS*(CI(1,i)^hCI))/(kCI^hCI+(CI(1,i)^hCI))+(betaHA*(CI(1,i)^hCI))/(kCI^hCI+(CI(1,i)^hCI));  %Creation of Salicylic acid
       
       ah= gammaAS*AS(1,i); %Degradation of Salicylic acid
       
       ai=eAS*AS(1,i);  %Export of Salicylic acid
       
       %-----------------------------------------------------------------%
       %                         Gillespie Algorithm
       %-----------------------------------------------------------------%
       
       med=a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad+ae+af+ag+ah+ai;
       
       T(1,i+1)=T(1,i)+1/(med) * log(1/rand());     %Time distribution
       
       n=rand();
       
       
       if (n>0) && (n<a/med)      %A molecule of C is created
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i)+1;
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
           
       elseif (n>a/med) && (n<(a+b)/med) && (C(1,i)>0)   %A molecule of C is consumed
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i)-1;
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b)/med) && (n<(a+b+c)/med) && C(1,i)>0 %The complex C-S is created
           
           S(1,i)=Stotal-1;
           C(1,i+1)=C(1,i)-1;
           CS(1,i+1)=CS(1,i)+1;
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c)/med) && (n<(a+b+c+d)/med) && (CS(1,i)>0) && QQ(1,i)>0 %The complex C-S attaches to chitin and frees the sensor
           
           S(1,i)=Stotal+1;
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i)-1;
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i)-1;
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d)/med) && (n<(a+b+c+d+e)/med) && CS(1,i)>0 && QQ(1,i)>0 %The complex CS is destroyed
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i)-1;
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e)/med) && (n<(a+b+c+d+e+f)/med)  %Chitoporin is created
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i)+1;
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f)/med) && (n<(a+b+c+d+e+f+g)/med) && P(1,i)>0 %The chitoporine is diluted
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i)-1;
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f+g)/med) && (n<(a+b+c+d+e+f+g+h)/med) % Thei chitinase is created
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i)+1;
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f+g+h)/med) && (n<(a+b+c+d+e+f+g+h+k)/med) && Ai(1,i)>0 %The chitinase inside the cell is destructed
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i)-1;
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f+g+h+k)/med) && (n<(a+b+c+d+e+f+g+h+k+l)/med) && (Ai(1,i)>0) % The chitinase is exported
           
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i)-1;
           Ao(1,i+1)=Ao(1,i)+1;
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f+g+h+k+l)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m)/med)&& Ao(1,i)>0 %the hitinase outside the cell is diluted
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i)-1;
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f+g+h+k+l+m)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o)/med) && (Ao(1,i)>0) && QQ(1,i)>0 %The chitinase clivajes the chitin
           
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i)-1;
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p)/med) && QQ(1,i)>0 %Thi chitin enters to cell
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i)+1;
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q)/med)  %The LuxI concentration inside the cell increases
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i)+1;
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
           
       elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r)/med) && Io(1,i)>0 %Import of LuxI
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i)+1;
           Io(1,i+1)=Io(1,i)-1;
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s)/med) && (Ii(1,i)>0)  %The LuxI concentration inside the cell decreases
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i)-1;
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t)/med) && Ii(1,i)>0 %Export of LuxI
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i)-1;
           Io(1,i+1)=Io(1,i)+1;
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u)/med) && (Ii(1,i)>0) && R(1,i)>0 %The complex LuxI-LuR is created
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i)-1;
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i)+1;
           R(1,i+1)=R(1,i)-1;
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v)/med) &&Io(1,i)>0 %Dilution of LuxI outsed the cell
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i)-1;
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w)/med) && (IR(1,i)>0)  %LuxI-LuxR is destryed
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i)-1;
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x)/med) % LuxR is created
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i)+1;
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y)/med) && (R(1,i)>0)  %R decresases
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i)-1;
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z)/med) %CI is created
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i)+1;
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa)/med) && (CI(1,i)>0)  %CI decreases
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i)-1;
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab)/med) %HB increases
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i)+1;
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac)/med) && (HB(1,i)>0)  %HB decreases
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i)-1;
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad)/med) && HA(1,i)>1 && HB(1,i)>1 %Inactivation of the toxin
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i)-2;
           HA(1,i+1)=HA(1,i)-2;
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad+ae)/med) && (HA(1,i)>0)  %HA decreases
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i)-1;
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad+ae)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad+ae+af)/med)   %HA decreases
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i)+1;
           AS(1,i+1)=AS(1,i);
           
       elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad+ae+af)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad+ae+af+ag)/med)   %AS Increases
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i)+1;
           
       elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad+ae+af+ag)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad+ae+af+ag+ah)/med) && AS(1,i)>0  %AS Decreases
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i)-1;
           
       elseif (n>(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad+ae+af+ag+ah)/med) && (n<(a+b+c+d+e+f+g+h+k+l+m+o+p+q+r+s+t+u+v+w+x+y+z+aa+ab+ac+ad+ae+af+ag+ah+ai)/med) && AS(1,i)>0 %AS Decreases
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i)-1;
           
       else
           
           S(1,i+1)= S(1,i);
           C(1,i+1)=C(1,i);
           CS(1,i+1)=CS(1,i);
           P(1,i+1)=P(1,i);
           Ai(1,i+1)=Ai(1,i);
           Ao(1,i+1)=Ao(1,i);
           Q(1,i+1)=Q(1,i);
           Ii(1,i+1)=Ii(1,i);
           Io(1,i+1)=Io(1,i);
           IR(1,i+1)=IR(1,i);
           R(1,i+1)=R(1,i);
           CI(1,i+1)=CI(1,i);
           HB(1,i+1)=HB(1,i);
           HA(1,i+1)=HA(1,i);
           AS(1,i+1)=AS(1,i);
           
           
           
           
           
       end
       
       if T(1,i) >= Tmaximo
           
           break
           
       end
       
       
   end
   
   %---------------------------------------------------------------------%
   %                       Converting to regular range
   %---------------------------------------------------------------------%
   
   tmax=Tmaximo;%max(T);
   stepp=0.001;
   
   [TFixed, CFixed]=regintervalfixed(T,C,stepp,tmax);
   [TFixed, CSFixed]=regintervalfixed(T,CS,stepp,tmax);
   [TFixed, PFixed]=regintervalfixed(T,P,stepp,tmax);
   [TFixed, AiFixed]=regintervalfixed(T,Ai,stepp,tmax);
   [TFixed, AoFixed]=regintervalfixed(T,Ao,stepp,tmax);
   [TFixed, QFixed]=regintervalfixed(T,Q,stepp,tmax);
   [TFixed, IiFixed]=regintervalfixed(T,Ii,stepp,tmax);
   [TFixed, IoFixed]=regintervalfixed(T,Io,stepp,tmax);
   [TFixed, IRFixed]=regintervalfixed(T,IR,stepp,tmax);
   [TFixed, RFixed]=regintervalfixed(T,R,stepp,tmax);
   [TFixed, CIFixed]=regintervalfixed(T,CI,stepp,tmax);
   [TFixed, HBFixed]=regintervalfixed(T,HB,stepp,tmax);
   [TFixed, HAFixed]=regintervalfixed(T,HA,stepp,tmax);
   [TFixed, ASFixed]=regintervalfixed(T,AS,stepp,tmax);
   [TFixed, QQFixed]=regintervalfixed(T,QQ,stepp,tmax);
   [TFixed, SFixed]=regintervalfixed(T,S,stepp,tmax);
   
   for k=1:length(TFixed)
       
       
       TMatrix(j,k)=TFixed(1,k);
       CMatrix(j,k)=CFixed(1,k);
       CSMatrix(j,k)=CSFixed(1,k);
       PMatrix(j,k)=PFixed(1,k);
       AiMatrix(j,k)=AiFixed(1,k);
       AoMatrix(j,k)=AoFixed(1,k);
       QMatrix(j,k)=QFixed(1,k);
       IiMatrix(j,k)=IiFixed(1,k);
       IoMatrix(j,k)=IoFixed(1,k);
       IRMatrix(j,k)=IRFixed(1,k);
       RMatrix(j,k)=RFixed(1,k);
       CIMatrix(j,k)=CIFixed(1,k);
       HBMatrix(j,k)=HBFixed(1,k);
       HAMatrix(j,k)=HAFixed(1,k);
       ASMatrix(j,k)=ASFixed(1,k);
       QQMatrix(j,k)=QQFixed(1,k);
       SMatrix(j,k)=SFixed(1,k);
       
       
   end
   
   
end
%
%-------------------------------------------------------------------------%
%                           Mean values for answers
%-------------------------------------------------------------------------%
%
%
%
%
%
TPromedio=mean(TMatrix);
CPromedio=mean(CMatrix);
CSPromedio=mean(CSMatrix);
PPromedio=mean(PMatrix);
AiPromedio=mean(AiMatrix);
AoPromedio=mean(AoMatrix);
QPromedio=mean(QMatrix);
IiPromedio=mean(IiMatrix);
IoPromedio=mean(IoMatrix);
IRPromedio=mean(IRMatrix);
RPromedio=mean(RMatrix);
CIPromedio=mean(CIMatrix);
HBPromedio=mean(HBMatrix);
HAPromedio=mean(HAMatrix);
ASPromedio=mean(ASMatrix);
QQPromedio=mean(QQMatrix);
SPromedio=mean(SMatrix);
%
%-------------------------------------------------------------------------%
%                       Standard deviation for answers
%-------------------------------------------------------------------------%
%
CDesv=std(CMatrix);
CSDesv=std(CSMatrix);
PDesv=std(PMatrix);
AiDesv=std(AiMatrix);
AoDesv=std(AoMatrix);
QDesv=std(QMatrix);
IiDesv=std(IiMatrix);
IoDesv=std(IoMatrix);
IRDesv=std(IRMatrix);
RDesv=std(RMatrix);
CIDesv=std(CIMatrix);
HBDesv=std(HBMatrix);
HADesv=std(HAMatrix);
ASDesv=std(ASMatrix);
QQDesv=std(QQMatrix);
SDesv=std(SMatrix);
%
%-------------------------------------------------------------------------%
%                                 Plotting
%-------------------------------------------------------------------------%
%
plot(TPromedio,ASPromedio,TPromedio,CIPromedio,TPromedio,QQPromedio/50);
xlabel('Time')
%
ylabel('Number of Molecules')
%legend('CBP','CBP-S','Chitoporin','Chitinase inside','Chitinase outside','Chitin monomers','Quitin');
%
%
%
%
%
%
%
%
Retrieved from "http://2012.igem.org/Team:Colombia/Modeling/Scripting"