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 % %