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

Sensitivity Analysis
This code points out which are the critical parameters in the system (those that change the response drastically).
% 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=5; % Dimerization rate of TetR ---[Sensitivity Analysis]--- [5-50/5] % 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=1; % TR basal production rate ---[Sensitivity Analysis]--- [0.5-5/0.5] ata=0.1; % TA basal production rate ---[Sensitivity Analysis]--- [0.1-10/0.5] % 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=5; % Maximum rate of TR expression ---[Sensitivity Analysis]--- [1-10/1] bta=5; % Maximum rate of TA expression ---[Sensitivity Analysis]--- [1-10/1] % ho=1.5; % LuxO*- DNA coupling rate ---[Sensitivity Analysis]--- [0.5-2.5/0.5] htr=5; % TRdomain-DNA coupling rate ---[Sensitivity Analysis]--- [1-10/1] % n=2; % Hill's coefficient ---[Sensitivity Analysis]--- [1-4/0.5] % h=400; % Maximum time m=1/10; % Step lenght [h] % param= (5:50:5); % Range for the parameter to vary, evaluate at least 2 orders of magnitude chose maximum 20 points. rta=zeros(3,length(param)); % Array that saves the canche in concentration of LuxOF, the response and time of response. for w=1:length(param) % Varying the value of the parameter kttr= param(w); % New parameter value (This MUST be changed everytime a parameter is being evaluated) % % %--------------------------------------------% % Runge Kutta 4 order aproximation % %--------------------------------------------% % % % 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(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 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 % R=x(:,10); Of=x(:,5);% % % it=141; % for it=1011:3001 % % if (R(it)-R(1011))/(R(1011))> 1.5 % If the response changes 1.5-fold relative to previous iteration % rta(1,w)= l(it)-100; % % % end % % % end % % % % rta(2,w)= Of(2701)-Of(801); rta(3,w)= R(2701)-R(801); % end % %Here Labels should be changed for them to make proper sense for each evaluated parameter % figure (1) plot(param,rta(1,:),'r') xlabel('Dimerization rate of TetR'); ylabel('Time of response') title('Sensitivity analysis') % figure (2) plot(param,rta(2,:),'g') xlabel('Dimerization rate of TetR'); ylabel('Change in LuxO Phosphorylated') title('Sensitivity analysis') % figure (3) plot(param,rta(3,:)) xlabel('Dimerization rate of TetR'); ylabel('Change in response (nM)') title('Sensitivity analysis')

Stochastic Model
Sometimes probabilistic models better describe certain systems