Team:Colombia/Scripting
From 2014.igem.org
(Difference between revisions)
Line 188: | Line 188: | ||
% | % | ||
% | % | ||
+ | </textarea> | ||
+ | <b><font color="#8A0808" size="3" >Metropolis-Hastings Algorithm</font></b> | ||
+ | <br> | ||
+ | This code determines which set of values for missing parameters yield the most desirable response. | ||
+ | <textarea> | ||
+ | function y=model(patr, pbtr, pho, pkttr, pata, pbta, phtr,pn) | ||
+ | global kcc kcd kcu kuc kuo kou kttr gcs gcsa guf gu gof go gtr gta gr gttr acs au ao ar atr ata bcu buc btr buo bou bta ho htr n | ||
+ | |||
+ | % | ||
+ | %--------------------------------------------% | ||
+ | % PARAMETERS % | ||
+ | %--------------------------------------------% | ||
+ | % | ||
+ | kcc=1; % CAI1 and CqsS coupling Rate | ||
+ | kcd=1; % CAI1 and CqsS decoupling Rate | ||
+ | kcu=2; % Michaelis Menten Constant for the phoshporylation CqsS-LuxU | ||
+ | kuc=1; % Michaelis Menten Constant of the phosphase process CqSS-LuxU | ||
+ | kuo=2; % Michaelis Menten Constant for the phoshporylation LuxU-LuxO | ||
+ | kou=1; % Michaelis Menten Constant of the phosphase process LuxU-LuxO | ||
+ | kttr=pkttr; %Dimerization rate of TetR | ||
+ | |||
+ | gcs=0.5; % CqsS protein decay rate | ||
+ | gcsa=1; % CqsS* protein decay rate | ||
+ | guf=0.12; % LuxU* protein decay rate | ||
+ | gu=0.65; % LuxU protein decay rate | ||
+ | gof=0.12; % LuxO* protein decay rate | ||
+ | go=0.12; % LuxO protein decay rate | ||
+ | gtr=0.12; % Ptet Repressor protein decay rate | ||
+ | gttr=0.12; %Ptet2 dimer Repressor decay rate | ||
+ | gta=0.12; % Ptet Activator protein decay rate | ||
+ | gr=0.12; % Response molecule decay rate | ||
+ | % | ||
+ | acs=5; % CS basal production rate | ||
+ | au=5; % LuxU basal production rate | ||
+ | ao=5; % LuxO basal production rate | ||
+ | ar=0; % response molecule basal production rate | ||
+ | atr=patr; % TR basal production rate | ||
+ | ata=pata; % TA basal production rate | ||
+ | % | ||
+ | bcu=1; % Phosphorylation rate of CqsS to LuxU | ||
+ | buc=0.1; % Phosphase rate of CqsS to LuxU | ||
+ | buo=3.2; % Phosphorylation rate LuxU-LuxO | ||
+ | bou=1.58; % Phosphase rate LuxU-LuxO | ||
+ | btr=pbtr; % Maximum rate of TR expression | ||
+ | bta=pbta; % Maximum rate of TA expression | ||
+ | % | ||
+ | ho=pho; % LuxO*- DNA coupling rate | ||
+ | htr=phtr; % TRdomain-DNA coupling rate | ||
+ | % | ||
+ | n=pn;% Hill coefficient | ||
+ | % | ||
+ | % | ||
+ | |||
+ | |||
+ | % | ||
+ | %--------------------------------------------% | ||
+ | % Runge Kutta 4 order aproximation % | ||
+ | %--------------------------------------------% | ||
+ | % | ||
+ | h=210; % Maximum time | ||
+ | m=1/12; % Step lenght [h] | ||
+ | conInd=[1,0,1,1,1,1,0,1,1,1]; % | ||
+ | l=(0:m:h)'; % Time vector | ||
+ | x=zeros(length(l),length(conInd)); % Variable's Matrix, the variables in each column and time in each row | ||
+ | C=zeros(1,length(l)); % zero vector to be filled | ||
+ | x(1,:)=conInd; % | ||
+ | % | ||
+ | for k=1:length(l)-1 | ||
+ | xk=x(k,:); % Capture of the last value in the matrix, the actual variable values | ||
+ | k1=difeq_Metropolis(l(k),xk); % First slope of the RK4 method | ||
+ | k2=difeq_Metropolis(l(k)+m/2,xk+(m/2*k1)'); % Second slope of the RK4 method | ||
+ | k3=difeq_Metropolis(l(k)+m/2,xk+(m/2*k2)'); % Third slope of the RK4 method | ||
+ | k4=difeq_Metropolis(l(k)+m,xk+(m*k3)'); % Fourth slope of the RK4 method | ||
+ | xk1=xk+m/6*(k1+2*k2+2*k3+k4)'; % New value calculation for the variables | ||
+ | %xk1=xk+m*difeq_Metropolis(l(k),xk)'; % Another aproximation, Newton's method | ||
+ | xk2=zeros(1,length(xk1)); % | ||
+ | % | ||
+ | for p=1:length(xk1) | ||
+ | if(xk1(p)<0.00000001) | ||
+ | xk2(p)=0; | ||
+ | else | ||
+ | xk2(p)=xk1(p); | ||
+ | end | ||
+ | end | ||
+ | x(k+1,:)=xk2; % Actualizacion del nuevo vector de variables en la matriz | ||
+ | end | ||
+ | % | ||
+ | for j=1:length(l) | ||
+ | if (l(j)>100 && l(j)<200) | ||
+ | C(j)=20; | ||
+ | else | ||
+ | C(j)=0; | ||
+ | |||
+ | end | ||
+ | end | ||
+ | % | ||
+ | |||
+ | CS=x(:,1); % | ||
+ | CSa=x(:,2); % | ||
+ | Uf=x(:,3); % | ||
+ | U=x(:,4); % | ||
+ | Of=x(:,5); % | ||
+ | O=x(:,6); % | ||
+ | TR=x(:,7); % | ||
+ | TTR=x(:,8); % | ||
+ | TA=x(:,9); % | ||
+ | R=x(:,10); % | ||
+ | |||
+ | %figure; % | ||
+ | %plot(l,C,l,R,'r'); % | ||
+ | %legend ('CAI','Chromoprotein'); % | ||
+ | %xlabel('Time [h]'); % | ||
+ | %ylabel('Concentration [nanomolar]'); % | ||
+ | %title('Response'); % | ||
+ | |||
+ | y=[R(95*12),R(100*12),R(105*12),R(110*12),R(115*12),R(120*12),R(125*12),R(130*12),R(135*12),R(140*12),R(145*12),R(150*12),R(155*12),R(160*12),R(165*12),R(170*12)]; | ||
+ | % | ||
+ | % | ||
+ | RFP=[3,3.1,3.15,3.5,4.8,6,11,15,16.5,18,18.5,19,19.5,20,20.25,20.5]; | ||
+ | t=[95,100,105,110,115,120,125,130,135,140,145,150,155,160,165,170]; | ||
+ | |||
+ | steps=5000; | ||
+ | |||
+ | atr_walk=zeros(1,steps+1); | ||
+ | btr_walk=zeros(1,steps+1); | ||
+ | ho_walk=zeros(1,steps+1); | ||
+ | kttr_walk=zeros(1,steps+1); | ||
+ | ata_walk=zeros(1,steps+1); | ||
+ | bta_walk=zeros(1,steps+1); | ||
+ | htr_walk=zeros(1,steps+1); | ||
+ | n_walk=zeros(1,steps+1); | ||
+ | |||
+ | atr_walk(1)=random('unif',0.5,1.5); | ||
+ | btr_walk(1)=random('unif',4,6); | ||
+ | ho_walk(1)=random('unif',1.2,1.8); | ||
+ | kttr_walk(1)=random('unif',2,7); | ||
+ | ata_walk(1)=random('unif',0.1,0.5); | ||
+ | bta_walk(1)=random('unif',4,7); | ||
+ | htr_walk(1)=random('unif',4,7); | ||
+ | n_walk(1)=random('unif',1,3); | ||
+ | |||
+ | for i=1:steps | ||
+ | |||
+ | atr_prime=random('unif',atr_walk(i)-0.2,atr_walk(i)+0.2); | ||
+ | btr_prime=random('unif',btr_walk(i)-0.2,btr_walk(i)+0.2); | ||
+ | ho_prime=random('unif',ho_walk(i)-0.2,ho_walk(i)+0.2); | ||
+ | kttr_prime=random('unif',kttr_walk(i)-0.2,kttr_walk(i)+0.2); | ||
+ | ata_prime=random('unif',ata_walk(i)-0.2,ata_walk(i)+0.2); | ||
+ | bta_prime=random('unif',bta_walk(i)-0.2,bta_walk(i)+0.2); | ||
+ | htr_prime=random('unif',htr_walk(i)-0.2,htr_walk(i)+0.2); | ||
+ | n_prime=random('unif',n_walk(i)-0.2,n_walk(i)+0.2); | ||
+ | |||
+ | model_i=model(atr_walk(i),btr_walk(i),ho_walk(i),kttr_walk(i),ata_walk(i),bta_walk(i),htr_walk(i),n_walk(i)); | ||
+ | model_prime=model(atr_prime,btr_prime,ho_prime,kttr_prime,ata_prime,bta_prime,htr_prime,n_prime); | ||
+ | |||
+ | alpha=likelihood_lineal(RFP,model_prime)/likelihood_lineal(RFP,model_i); | ||
+ | |||
+ | if(alpha>=1) | ||
+ | atr_walk(i+1)=atr_prime; | ||
+ | btr_walk(i+1)=btr_prime; | ||
+ | ho_walk(i+1)=ho_prime; | ||
+ | kttr_walk(i+1)=kttr_prime; | ||
+ | ata_walk(i+1)=ata_prime; | ||
+ | bta_walk(i+1)=bta_prime; | ||
+ | htr_walk(i+1)=htr_prime; | ||
+ | n_walk(i+1)=n_prime; | ||
+ | else | ||
+ | % beta=random('unif',0,1); | ||
+ | %if(beta<=alpha) | ||
+ | % atr_walk(i+1)=atr_prime; | ||
+ | % btr_walk(i+1)=btr_prime; | ||
+ | %ho_walk(i+1)=ho_prime; | ||
+ | %kttr_walk(i+1)=kttr_prime; | ||
+ | %ata_walk(i+1)=ata_prime; | ||
+ | %bta_walk(i+1)=bta_prime; | ||
+ | %htr_walk(i+1)=htr_prime; | ||
+ | %n_walk(i+1)=n_prime; | ||
+ | %else | ||
+ | atr_walk(i+1)=atr_walk(i); | ||
+ | btr_walk(i+1)=btr_walk(i); | ||
+ | ho_walk(i+1)=ho_walk(i); | ||
+ | kttr_walk(i+1)=kttr_walk(i); | ||
+ | ata_walk(i+1)=ata_walk(i); | ||
+ | bta_walk(i+1)=bta_walk(i); | ||
+ | htr_walk(i+1)=htr_walk(i); | ||
+ | n_walk(i+1)=n_walk(i); | ||
+ | %end | ||
+ | end | ||
+ | i | ||
+ | end | ||
+ | |||
+ | figure (1) | ||
+ | plot(t,RFP,'b',t,model(atr_walk(steps+1),btr_walk(steps+1),ho_walk(steps+1),kttr_walk(steps+1),ata_walk(steps+1),bta_walk(steps+1),htr_walk(steps+1),n_walk(steps+1)),'r') | ||
+ | xlabel('tiempo') | ||
+ | ylabel('Reportero') | ||
+ | title('Respuesta') | ||
+ | |||
+ | figure (2) | ||
+ | scatter((1:steps+1),atr_walk) | ||
+ | xlabel('step') | ||
+ | ylabel('valor') | ||
+ | title('atr') | ||
+ | |||
+ | figure (3) | ||
+ | scatter((1:steps+1),btr_walk) | ||
+ | xlabel('step') | ||
+ | ylabel('valor') | ||
+ | title('btr') | ||
+ | |||
+ | figure (4) | ||
+ | scatter((1:steps+1),ho_walk) | ||
+ | xlabel('step') | ||
+ | ylabel('valor') | ||
+ | title('ho') | ||
+ | |||
+ | figure (5) | ||
+ | scatter((1:steps+1),kttr_walk) | ||
+ | xlabel('step') | ||
+ | ylabel('valor') | ||
+ | title('kttr') | ||
+ | |||
+ | |||
+ | figure (6) | ||
+ | scatter((1:steps+1),ata_walk) | ||
+ | xlabel('step') | ||
+ | ylabel('valor') | ||
+ | title('ata') | ||
+ | |||
+ | |||
+ | figure (7) | ||
+ | scatter((1:steps+1),bta_walk) | ||
+ | xlabel('step') | ||
+ | ylabel('valor') | ||
+ | title('bta') | ||
+ | |||
+ | |||
+ | figure (8) | ||
+ | scatter((1:steps+1),htr_walk) | ||
+ | xlabel('step') | ||
+ | ylabel('valor') | ||
+ | title('htr') | ||
+ | |||
+ | |||
+ | figure (9) | ||
+ | scatter((1:steps+1),n_walk) | ||
+ | xlabel('step') | ||
+ | ylabel('valor') | ||
+ | title('n') | ||
+ | % | ||
+ | % | ||
+ | function y=input_Metropolis(t) | ||
+ | %--------------------------------------------% | ||
+ | % CAI-1 (C) pulse simulation % | ||
+ | %--------------------------------------------% | ||
+ | if t>100 && t<200 | ||
+ | y=20; | ||
+ | else | ||
+ | y=0; | ||
+ | end | ||
</textarea> | </textarea> | ||
</p> | </p> | ||
</div> | </div> | ||
</html> | </html> |
Revision as of 03:19, 16 October 2014
Scripting
Feel free to expand and scroll through the text boxes in order to further examine the code.
DETERMINISTIC MODEL
This code creates the differential equations governing the concentration dinamics of each protein in our model, finds the steady state solutions and then solves them using the numerical aproximation method Runge-Kutta Metropolis-Hastings Algorithm
This code determines which set of values for missing parameters yield the most desirable response.