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
%Stochastic Simulation for the system % %%----Parameters of the System---- global kcc kcd kcu kuc kuo kou kttr gcs gcsa guf gu gof go gtr gta gr acs au ao ar atr ata bcu buc btr buo bou bta ho htr n gttr % % numCells = 20; numEvents= 50000; % % % % kcc=1; % CAI1 and CqsS coupling Rate kcd=1; % CAI1 and CqsS decoupling Rate kcu=1.2044000; % Michaelis Menten Constant for the phoshporylation CqsS-LuxU kuc=0.6022000; % Michaelis Menten Constant of the phosphase process CqSS-LuxU kuo=1.2044000; % Michaelis Menten Constant for the phoshporylation LuxU-LuxO kou=0.6022000; % Michaelis Menten Constant of the phosphase process LuxU-LuxO kttr=5; %Dimerization rate of TetR % gcs=0.500; % CqsS protein decay rate gcsa=1; % CqsS* protein decay rate guf=0.120; % LuxU* protein decay rate gu=0.650; % LuxU protein decay rate gof=0.120; % LuxO* protein decay rate go=0.650; % LuxO protein decay rate gtr=0.120; % Ptet Repressor protein decay rate gttr=0.120; %Ptet2 dimer Repressor decay rate gta=0.120; % Ptet Activator protein decay rate gr=0.050; % Response molecule decay rate % acs=5; % CS basal production rate au=5; % LuxU basal production rate ao=5; % LuxO basal production rate ar=0.1; % response molecule basal production rate atr=1; % TR basal production rate ata=0.1; % 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=5; % Maximum rate of TR expression bta=5; % Maximum rate of TA expression % ho=1.5; % LuxO*- DNA coupling rate htr=5; % TRdomain-DNA coupling rate % n=3;% Hill coefficient % %%----Stochastic Simulation---- % %Variables of the System C=zeros(numCells,numEvents); CS=zeros(numCells,numEvents); CSa=zeros(numCells,numEvents); Uf=zeros(numCells,numEvents); U=zeros(numCells,numEvents); Of=zeros(numCells,numEvents); O=zeros(numCells,numEvents); TR=zeros(numCells,numEvents); TRR=zeros(numCells,numEvents); TA=zeros(numCells,numEvents); R=zeros(numCells,numEvents); % C(:,1)=0; CS(:,1)=10; CSa(:,1)=1; Uf(:,1)=15; U(:,1)=5; Of(:,1)=40; O(:,1)=5; TA(:,1)=10; TR(:,1)=2; TRR(:,1)=1; % % % % t = zeros(numCells,numEvents); %Time tr = zeros(numCells,numEvents); % % %Stochastic Simulation % for i =1:numCells i for j = 1:numEvents j %Random numbers drawn from uniform distribution between zero and %one rand1 = rand(1); %This determine time rand2 = rand(1); %This determine Event %Events (Description in attached document) e1= kcc*CS(i,j)*C(i,j); e2= kcd*CSa(i,j); e3= gcs*CS(i,j); e4= acs; e5= gcsa*CSa(i,j); e6= au; e7= (bcu*U(i,j)*CS(i,j))/(kcu+U(i,j)); e8= (buc*Uf(i,j)*CSa(i,j))/(kuc+Uf(i,j)); e9= (buo*O(i,j)*Uf(i,j))/(kuo+O(i,j)); e10=(bou*Of(i,j)*U(i,j))/(kou+Of(i,j)); e11=gu*U(i,j); e12=guf*Uf(i,j); e13=ao; e14=go*O(i,j); e15=gof*Of(i,j); e16=atr; e17=(btr*Of(i,j)^n)/((ho.^n)+(Of(i,j).^n)); e18=kttr*TR(i,j)^2; e19=gtr*TR(i,j); e20=ata; e21=ar; e22=bta/(1+(TRR(i,j)/TA(i,j))+(htr/TA(i,j))); e23=gta*TA(i,j); e24=gr*R(i,j); e25=TRR(i,j)*gttr; %Vector of Events vecEvents = [e1,e2,e3,e4,e5,e6,e7,e8,e9,e10,e11,e12,e13,e14,e15,e16,e17,e18,e19,e20,e21,e22,e23,e24,e25]; %Sum of Events sumEvents = sum(vecEvents); %Waiting time until the next ocurrence of a reaction time %Probabilistic Distribution tau= (1/sumEvents)*log(1/rand1); tr(i,j)=tau; t(i,j+1)= t(i,j)+ tau; %CAI outside the cell C(i,j+1)=funcImpulso(t(i,j+1)); %if(t(i,j)<130) % C(i,j)=funcImpulso(t(i,j+1))*957/2575; %end %Simulation Creation/Destruction of each molecule if (rand2<(vecEvents(1)/sumEvents)) %CqsS activation if (CS(i,j)<=0) CS(i,j+1)=0; else CS(i,j+1)= CS(i,j)-1; end CSa(i,j+1)= CSa(i,j)+1; U(i,j+1)= U(i,j); Uf(i,j+1)= Uf(i,j); O(i,j+1)= O(i,j); Of(i,j+1)= Of(i,j); TR(i,j+1)= TR(i,j); TRR(i,j+1)= TRR(i,j); TA(i,j+1)= TA(i,j); R(i,j+1)= R(i,j); elseif (rand2< (sum(vecEvents(1:2)))/sumEvents) %CqsS desactivation if (CSa(i,j)<=0) CSa(i,j+1)= 0; else CSa(i,j+1)= CSa(i,j)-1; end CS(i,j+1)= CS(i,j)+1; U(i,j+1)= U(i,j); Uf(i,j+1)= Uf(i,j); O(i,j+1)= O(i,j); Of(i,j+1)= Of(i,j); TR(i,j+1)= TR(i,j); TRR(i,j+1)= TRR(i,j); TA(i,j+1)= TA(i,j); R(i,j+1)= R(i,j); elseif (rand2< (sum(vecEvents(1:3)))/sumEvents) %CqsS degradation if (CS(i,j)<=0) CS(i,j+1)=0; else CS(i,j+1)= CS(i,j)-1; end CSa(i,j+1)= CSa(i,j); U(i,j+1)= U(i,j); Uf(i,j+1)= Uf(i,j); O(i,j+1)= O(i,j); Of(i,j+1)= Of(i,j); TR(i,j+1)= TR(i,j); TRR(i,j+1)= TRR(i,j); TA(i,j+1)= TA(i,j); R(i,j+1)= R(i,j); elseif (rand2< (sum(vecEvents(1:4)))/sumEvents) %CqsS basal production CS(i,j+1)= CS(i,j)+1; CSa(i,j+1)= CSa(i,j); U(i,j+1)= U(i,j); Uf(i,j+1)= Uf(i,j); O(i,j+1)= O(i,j); Of(i,j+1)= Of(i,j); TR(i,j+1)= TR(i,j); TRR(i,j+1)= TRR(i,j); TA(i,j+1)= TA(i,j); R(i,j+1)= R(i,j); elseif (rand2< (sum(vecEvents(1:5)))/sumEvents) %CqsS* degradation if (CSa(i,j)<=0) CSa(i,j+1)=0; else CSa(i,j+1)= CSa(i,j)-1; end CS(i,j+1)= CS(i,j); U(i,j+1)= U(i,j); Uf(i,j+1)= Uf(i,j); O(i,j+1)= O(i,j); Of(i,j+1)= Of(i,j); TR(i,j+1)= TR(i,j); TRR(i,j+1)= TRR(i,j); TA(i,j+1)= TA(i,j); R(i,j+1)= R(i,j); elseif (rand2< (sum(vecEvents(1:6)))/sumEvents) %LuxU basal production U(i,j+1)= U(i,j)+1; CS(i,j+1)= CS(i,j); CSa(i,j+1)= CSa(i,j); Uf(i,j+1)= Uf(i,j); O(i,j+1)= O(i,j); Of(i,j+1)= Of(i,j); TR(i,j+1)= TR(i,j); TRR(i,j+1)= TRR(i,j); TA(i,j+1)= TA(i,j); R(i,j+1)= R(i,j); elseif (rand2< (sum(vecEvents(1:7)))/sumEvents) %LuxU fosforilation if (U(i,j)<=0) U(i,j+1)=0; else U(i, j+1)= U(i,j)-1; end CS(i,j+1)= CS(i,j); CSa(i,j+1)= CSa(i,j); Uf(i,j+1)= Uf(i,j)+1; O(i,j+1)= O(i,j); Of(i,j+1)= Of(i,j); TR(i,j+1)= TR(i,j); TRR(i,j+1)= TRR(i,j); TA(i,j+1)= TA(i,j); R(i,j+1)= R(i,j); elseif (rand2< (sum(vecEvents(1:8)))/sumEvents) %LuxU dephosphorilation if (Uf(i,j)<=0) Uf(i,j+1)=0; else Uf(i, j+1)= Uf(i,j)-1; end CS(i,j+1)= CS(i,j); CSa(i,j+1)= CSa(i,j); U(i,j+1)= U(i,j)+1; O(i,j+1)= O(i,j); Of(i,j+1)= Of(i,j); TR(i,j+1)= TR(i,j); TRR(i,j+1)= TRR(i,j); TA(i,j+1)= TA(i,j); R(i,j+1)= R(i,j); elseif (rand2< (sum(vecEvents(1:9)))/sumEvents) %LuxO phosphorilation if (Uf(i,j)<=0 || O(i,j)<=0 ) if(Uf(i,j)<0) Uf(i,j+1)=0; else Uf(i,j+1)=Uf(i,j); end if(O(i,j)<0) O(i,j+1)=0; else O(i,j+1)=O(i,j); end else Uf(i, j+1)= Uf(i,j)-1; O(i,j+1)=O(i,j)-1; end CS(i,j+1)= CS(i,j); CSa(i,j+1)= CSa(i,j); U(i,j+1)= U(i,j)+1; Of(i,j+1)= Of(i,j)+1; TR(i,j+1)= TR(i,j); TRR(i,j+1)= TRR(i,j); TA(i,j+1)= TA(i,j); R(i,j+1)= R(i,j); elseif (rand2< (sum(vecEvents(1:10)))/sumEvents) %LuxO dephosphorilation if (U(i,j)<=0) U(i,j+1)=0; else U(i,j+1)= U(i,j)-1; end if (Of(i,j)<=0) Of(i,j+1)=0; else Of(i, j+1)= Of(i,j)-1; end CS(i,j+1)= CS(i,j); CSa(i,j+1)= CSa(i,j); Uf(i,j+1)= Uf(i,j)+1; O(i,j+1)= O(i,j)+1; TR(i,j+1)= TR(i,j); TRR(i,j+1)= TRR(i,j); TA(i,j+1)= TA(i,j); R(i,j+1)= R(i,j); elseif (rand2< (sum(vecEvents(1:11)))/sumEvents) %LuxU degradation if (U(i,j)<=0) U(i,j+1)=0; else U(i, j+1)= U(i,j)-1; end CS(i,j+1)= CS(i,j); CSa(i,j+1)= CSa(i,j); Uf(i,j+1)= Uf(i,j); O(i,j+1)= O(i,j); Of(i,j+1)= Of(i,j); TR(i,j+1)= TR(i,j); TRR(i,j+1)= TRR(i,j); TA(i,j+1)= TA(i,j); R(i,j+1)= R(i,j); elseif (rand2< (sum(vecEvents(1:12)))/sumEvents ) %LuxU* degradation if (Uf(i,j)<=0) Uf(i,j+1)=0; else Uf(i, j+1)= Uf(i,j)-1; end CS(i,j+1)= CS(i,j); CSa(i,j+1)= CSa(i,j); U(i,j+1)= U(i,j); O(i,j+1)= O(i,j); Of(i,j+1)= Of(i,j); TR(i,j+1)= TR(i,j); TRR(i,j+1)= TRR(i,j); TA(i,j+1)= TA(i,j); R(i,j+1)= R(i,j); elseif (rand2< (sum(vecEvents(1:13)))/sumEvents )%Basal production LuxO U(i, j+1)= U(i,j); CS(i,j+1)= CS(i,j); CSa(i,j+1)= CSa(i,j); Uf(i,j+1)= Uf(i,j); O(i,j+1)= O(i,j)+1; Of(i,j+1)= Of(i,j); TR(i,j+1)= TR(i,j); TRR(i,j+1)= TRR(i,j); TA(i,j+1)= TA(i,j); R(i,j+1)= R(i,j); elseif (rand2< (sum(vecEvents(1:14)))/sumEvents) %LuxO degradation if (O(i,j)<=0) O(i,j+1)=0; else O(i, j+1)= O(i,j)-1; end CS(i,j+1)= CS(i,j); CSa(i,j+1)= CSa(i,j); U(i,j+1)= U(i,j); Uf(i,j+1)= Uf(i,j); Of(i,j+1)= Of(i,j); TR(i,j+1)= TR(i,j); TRR(i,j+1)= TRR(i,j); TA(i,j+1)= TA(i,j); R(i,j+1)= R(i,j); elseif (rand2< (sum(vecEvents(1:15)))/sumEvents ) %LuxO* degradation if (Of(i,j)<=0) Of(i,j+1)=0; else Of(i, j+1)= Of(i,j)-1; end CS(i,j+1)= CS(i,j); CSa(i,j+1)= CSa(i,j); U(i,j+1)= U(i,j); O(i,j+1)= O(i,j); Uf(i,j+1)= Uf(i,j); TR(i,j+1)= TR(i,j); TRR(i,j+1)= TRR(i,j); TA(i,j+1)= TA(i,j); R(i,j+1)= R(i,j); elseif (rand2< (sum(vecEvents(1:16)))/sumEvents) %Basal Production tetR TR(i, j+1)= TR(i,j)+1; CS(i,j+1)= CS(i,j); CSa(i,j+1)= CSa(i,j); U(i,j+1)= U(i,j); O(i,j+1)= O(i,j); Of(i,j+1)= Of(i,j); Uf(i,j+1)= Uf(i,j); TRR(i,j+1)= TRR(i,j); TA(i,j+1)= TA(i,j); R(i,j+1)= R(i,j); elseif (rand2< (sum(vecEvents(1:17)))/sumEvents )%Active production tetR TR(i, j+1)= TR(i,j)+1; CS(i,j+1)= CS(i,j); CSa(i,j+1)= CSa(i,j); U(i,j+1)= U(i,j); O(i,j+1)= O(i,j); Of(i,j+1)= Of(i,j); Uf(i,j+1)= Uf(i,j); TRR(i,j+1)= TRR(i,j); TA(i,j+1)= TA(i,j); R(i,j+1)= R(i,j); elseif (rand2< (sum(vecEvents(1:18)))/sumEvents) %tetR dimerization if (TR(i,j)<=0) TR(i,j+1)=0; TRR(i, j+1)= TRR(i,j); elseif (TR(i,j)==1) TR(i, j+1)=1; TRR(i, j+1)= TRR(i,j); else TR(i, j+1)= TR(i,j)-2; TRR(i, j+1)= TRR(i,j)+1; end CS(i,j+1)= CS(i,j); CSa(i,j+1)= CSa(i,j); U(i,j+1)= U(i,j); O(i,j+1)= O(i,j); Of(i,j+1)= Of(i,j); Uf(i,j+1)= Uf(i,j); TA(i,j+1)= TA(i,j); R(i,j+1)= R(i,j); elseif (rand2< (sum(vecEvents(1:19)))/sumEvents) %tetR degradation if (TR(i,j)<=0) TR(i,j+1)=0; else TR(i,j+1)= TR(i,j)-1; end CS(i,j+1)= CS(i,j); CSa(i,j+1)= CSa(i,j); U(i,j+1)= U(i,j); O(i,j+1)= O(i,j); Of(i,j+1)= Of(i,j); Uf(i,j+1)= Uf(i,j); TRR(i,j+1)= TRR(i,j); TA(i,j+1)= TA(i,j); R(i,j+1)= R(i,j); elseif (rand2< (sum(vecEvents(1:20)))/sumEvents )%Basal production Tet activator TA(i, j+1)= TA(i,j)+1; CS(i,j+1)= CS(i,j); CSa(i,j+1)= CSa(i,j); U(i,j+1)= U(i,j); O(i,j+1)= O(i,j); Of(i,j+1)= Of(i,j); Uf(i,j+1)= Uf(i,j); TRR(i,j+1)= TRR(i,j); TR(i,j+1)= TR(i,j); R(i,j+1)= R(i,j); elseif (rand2< (sum(vecEvents(1:21)))/sumEvents) %Basal production Reporter protein R(i, j+1)= R(i,j)+1; CS(i,j+1)= CS(i,j); CSa(i,j+1)= CSa(i,j); U(i,j+1)= U(i,j); O(i,j+1)= O(i,j); Of(i,j+1)= Of(i,j); Uf(i,j+1)= Uf(i,j); TRR(i,j+1)= TRR(i,j); TR(i,j+1)= TR(i,j); TA(i,j+1)= TA(i,j); elseif (rand2< (sum(vecEvents(1:22)))/sumEvents )%Active production Reporter protein & activator TA(i,j+1)= TA(i,j)+1; R(i,j+1)= R(i,j)+1; CS(i,j+1)= CS(i,j); CSa(i,j+1)= CSa(i,j); U(i,j+1)= U(i,j); O(i,j+1)= O(i,j); Of(i,j+1)= Of(i,j); Uf(i,j+1)= Uf(i,j); TRR(i,j+1)= TRR(i,j); TR(i,j+1)= TR(i,j); elseif (rand2< (sum(vecEvents(1:23)))/sumEvents) %Tet activator degradation if (TA(i,j)<=0) TA(i,j+1)=0; else TA(i, j+1)= TA(i,j)-1; end CS(i,j+1)= CS(i,j); CSa(i,j+1)= CSa(i,j); U(i,j+1)= U(i,j); O(i,j+1)= O(i,j); Of(i,j+1)= Of(i,j); Uf(i,j+1)= Uf(i,j); TRR(i,j+1)= TRR(i,j); TR(i,j+1)= TR(i,j); R(i,j+1)= R(i,j); elseif (rand2 < (sum(vecEvents(1:24)))/sumEvents) %Reporter protein activator degradation if (R(i,j)<=0) R(i,j+1)=0; else R(i, j+1)= R(i,j)-1; end CS(i,j+1)= CS(i,j); CSa(i,j+1)= CSa(i,j); U(i,j+1)= U(i,j); O(i,j+1)= O(i,j); Of(i,j+1)= Of(i,j); Uf(i,j+1)= Uf(i,j); TRR(i,j+1)= TRR(i,j); TR(i,j+1)= TR(i,j); TA(i,j+1)= TA(i,j); elseif (rand2 < (sum(vecEvents(1:25)))/sumEvents) %TetR dimers degradation if (TRR(i,j)<=0) TRR(i,j+1)=0; else TRR(i, j+1)= TRR(i,j)-1; end CS(i,j+1)= CS(i,j); CSa(i,j+1)= CSa(i,j); U(i,j+1)= U(i,j); O(i,j+1)= O(i,j); Of(i,j+1)= Of(i,j); Uf(i,j+1)= Uf(i,j); TR(i,j+1)= TR(i,j); TA(i,j+1)= TA(i,j); R(i, j+1)= R(i,j); end end end % Pro_CS=zeros(numEvents); Pro_CSa=zeros(numEvents); Pro_Uf=zeros(numEvents); Pro_U=zeros(numEvents); Pro_Of=zeros(numEvents); Pro_O=zeros(numEvents); Pro_TR=zeros(numEvents); Pro_TRR=zeros(numEvents); Pro_TA=zeros(numEvents); Pro_R=zeros(numEvents); % % for i=1:numEvents Pro_CS(i)=(sum(CS(:,i)))/numCells; end for i=1:numEvents Pro_CSa(i)=(sum(CSa(:,i)))/numCells; end for i=1:numEvents Pro_Uf(i)=(sum(Uf(:,i)))/numCells; end for i=1:numEvents Pro_U(i)=(sum(U(:,i)))/numCells; end for i=1:numEvents Pro_Of(i)=(sum(Of(:,i)))/numCells; end for i=1:numEvents Pro_O(i)=(sum(O(:,i)))/numCells; end for i=1:numEvents Pro_TR(i)=(sum(TR(:,i)))/numCells; end for i=1:numEvents Pro_TRR(i)=(sum(TRR(:,i)))/numCells; end for i=1:numEvents Pro_TA(i)=(sum(TA(:,i)))/numCells; end for i=1:numEvents Pro_R(i)=(sum(R(:,i)))/numCells; end % % % function [newTime newData] = regIntervalFixed2(Time,Data,step,timeMax) % % This function takes the output of a gillespie's simulation and turns it % into a matrix with a time vector of regular intervals. % % INPUT % Time: A matrix sized number of Simulations(Cells) by number of Events % that contains the time asociated with each event in each simulation. % Data: A matrix with the data. (events in columns and simulation in rows) % step: a scalar of the time interval. It must have the same units as the % entrances of the matrix Time. % timeMax: the maximum time required. It must have the same units as the % entrances of the matrix Time. % % OUTPUT: % newTime: a vector of time regularly spaced. % newData: a matrix where each row represents a simulation (cell) and each % column represents the number of molecules for each time in newTime L=size(Time,2); newTime=0:step:timeMax; newData = zeros(size(Data,1),length(newTime)); for k=1:size(Data,1) newdata=zeros(size(newTime)); j=1; for i=1:length(newTime) if step*(i-1)>Time(k,j) && j<L while step*(i-1)>Time(k,j) && j<L j=j+1; end end newdata(i)=Data(k,j); end newData(k,:) = newdata; end %newdata(L+1)=data(L); % % function impulse= funcImpulso(t) if (t>150 && t<310) impulse=100; % Molecules/L else impulse=0; %Molecules/L end % end % % function w=condiciones() % xi=[0.5 1000 5 1000]; % y=fsolve(@CondInd,xi,optimset('algorithm','levenberg-marquardt','maxiter',100000,'tolfun',1e-9)); % w=y; % end % %