"
Team:Colombia/Modeling/Scripting
From 2012.igem.org
(Difference between revisions)
(→RUST DIFFERENTIAL EQUATION SOLUTION) |
(→OPTIMIZATION) |
||
| Line 921: | Line 921: | ||
</div> | </div> | ||
| + | |||
| + | |||
| + | ==STOCHASTIC MODEL == | ||
Revision as of 03:28, 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;
