"
Team:Colombia/Modeling/Scripting
From 2012.igem.org
(Difference between revisions)
(→RUST DIFFERENTIAL EQUATION SOLUTION) |
(→RUST DIFFERENTIAL EQUATION SOLUTION) |
||
| Line 555: | Line 555: | ||
HA=x(:,13); | HA=x(:,13); | ||
AS=x(:,14); | AS=x(:,14); | ||
| - | + | % | |
% | % | ||
figure(1) | figure(1) | ||
Revision as of 13:35, 26 September 2012
Template:Https://2012.igem.org/User:Tabima
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')
%
