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 function y=det() % %--------------------------------------------% % Runge Kutta 4 order aproximation % %--------------------------------------------% % h=1200; % Maximum time m=0.005; % Step lenght [s] xi=[1,1,1,1,1,1,1,1,1]; % y=fsolve(@CondIni,xi); % conInd=[1,1,1,1,1,1,1,1,0]; % 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(l(k),xk); % First slope of the RK4 method k2=difeq(l(k)+m/2,xk+(m/2*k1)'); % Second slope of the RK4 method k3=difeq(l(k)+m/2,xk+(m/2*k2)'); % Third slope of the RK4 method k4=difeq(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(l(k),xk)'; % Another aproximation, Newton's method (comment line) 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; % Variable vector actualization in the matrix end % for j=1:length(l) if (l(j)>(300) & l(j)<(1000)) C(j)=2.25; 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); % TA=x(:,8); % R=x(:,9); % % % figure(1) % hold on; plot(l,R,'r') plot(l,C) plot(l,Of,'g') %plot(l,TR) hold off; xlabel('Time') ylabel('Concetration (micromolar)') title('Activator response') % end % % function y=difeq(t,x) % global kcc kcd dsu duo fsu fuo gcs gcsa guf gu gof go gtr gta gr acs au ao ar atr ata btr bta br ho htr n % %--------------------------------------------% % VARIABLES % %--------------------------------------------% % C=input(t); % Extracellular concentration of the cholerae autoinducer-1 (CAI-1) % CS=x(1); % Concentration of the membrane bound CqsS protein (CAI-1 unbound=Inactive) CSa=x(2); % Concentration of the membrane bound CqsS protein (CAI-1 bound=Active) Uf= x(3); % Phosphorelay protein LuxU (phosphorylated) Concentration U=x(4); % Phosphorelay protein LuxU (unphosphorylated) Concentration Of=x(5); % Phosphorelay protein LuxO (phosphorylated) Concentration O=x(6); % Phosphorelay protein LuxO (unphosphorylated) Concentration TR=x(7); % Ptet Repressor protein concentration TA=x(8); % Ptet Activator protein concentration R=x(9); % Response molecule concentration % %--------------------------------------------% % PARAMETERS % %--------------------------------------------% % kcc=1; % CAI1 and CqsS coupling Rate kcd=1; % CAI1 and CqsS decoupling Rate dsu=4; % LuxU* dephosphorylation rate through CqsS* duo=4; % LuxO* dephosphorylation rate through LuxU fsu=2; % LuxU* phosphorylation rate through CqsS fuo=2; % LuxO* phosphorylation rate through LuxU* % gcs=1; % CqsS protein decay rate gcsa=1; % CqsS* protein decay rate guf=1; % LuxU* protein decay rate gu=1; % LuxU protein decay rate gof=1; % LuxO* protein decay rate go=1; % LuxO protein decay rate gtr=1; % Ptet Repressor protein decay rate gta=1; % Ptet Activator protein decay rate gr=1; % Response molecule decay rate % acs=3; % CS basal production rate au=3; % LuxU basal production rate ao=3; % LuxO basal production rate ar=0.01; % response molecule basal production rate atr=0.01; % TR basal production rate ata=0.01; % TA basal production rate btr=5; % Maximum rate of TR expression bta=5; % Maximum rate of TA expression br=5; % Maximum rate of response molecule expression ho=1.5; % LuxO*- DNA coupling rate htr=2; % TRdomain-DNA coupling rate % n=1; % Hill coefficient % %--------------------------------------------% % Equations % %--------------------------------------------% % dCS = acs + kcd*CSa - C*CS*kcc - gcs*CS; % Differential equation governing the change in the concentration of the membrane bound CqsS protein (CAI-1 unbound=Inactive)through time dCSa = -kcd*CSa + C*CS*kcc - gcsa*CSa; % Differential equation governing the change in the concentration of the membrane bound CqsS protein (CAI-1 bound=Active) through time dUf = U*Of*duo-Uf*O*fuo-CSa*Uf*dsu + CS*U*fsu - guf*Uf; % Differential equation governing the change in the phosphorelay protein LuxU (phosphorylated) Concentration through time dU = au-U*Of*duo+Uf*O*fuo - CS*U*fsu + CSa*Uf*dsu -gu*U; % Differential equation governing the change in the phosphorelay protein LuxU (unphosphorylated) Concentration through time dOf = Uf*O*fuo - U*Of*duo - gof*Of; % Differential equation governing the change in the phosphorelay protein LuxO (phosphorylated) Concentration through time dO = ao - Uf*O*fuo + U*Of*duo - go*O; % Differential equation governing the change in the phosphorelay protein LuxO (unphosphorylated) Concentration through time dTR = atr + (btr*Of^n)/(ho^n+Of^n) - gtr*TR; % Differential equation governing the change in the ptet Repressor protein concentration through time dTA = ata + bta/((1+TR/TA)+(htr/TA)) - gta*TA; % Differential equation governing the change in the ptet Activator protein concentration through time dR = ar + br/((1+TR/TA)+(htr/TA)) - gr*R; % Differential equation governing the change in the response molecule concentration through time % y(1)=dCS; % y(2)=dCSa; % y(3)=dUf; % y(4)=dU; % y(5)=dOf; % y(6)=dO; % y(7)=dTR; % y(8)=dTA; % y(9)=dR; % % y=y'; % Transposing the vector because of matlab language restrictions % end % % function y=input(t) %--------------------------------------------% % CAI-1 (C) pulse simulation % %--------------------------------------------% if t>50 && t<100 y=10; elseif t>250 && t<1500 y=20; else y=0; end % % Metropolis-Hastings Algorithm
This code determines which set of values for missing parameters yield the most desirable response. 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