Team:Colombia/Modeling/Scripting

From 2012.igem.org

Revision as of 22:30, 26 September 2012 by Cesarquintana c (Talk | contribs)

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 THE SYSTEM FOR %RASTONIA

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;