% 6.581 PS4 % Problem 2 % lambda.m function [dxdt, dxdx] = lambda(x) % % [dxdt, dxdx] = lambda(x) implemets a simplified model of the % the bacteriophage lamda lysis lysogeny switch we saw in class. % The input is the current state vector: % % x = [MCI; %CI mRNA % CI; %CI Monomer % CI2; %CI Dimer % CI2B; %CI Dimer Bound to DNA % MCRO; %CRO mRNA % CRO; %CRO Monomer % CRO2; %CRO Dimer % CRO2B %CRO Dimer bound to DNA % ] % % The return is the the time derivative of x and the jacobian of x % % You will have to fill in the missing lines of code labled with % %% CHANGE ME % % One hint for finding fixed points: % Try starting the system with CRO=1 but CI=0. Where does the system go? % Try starting the system with CRO=0 but CI=1. Where does the system go? % Is the system stable at the end of the run? % % There is one other fixed point? Do you expect this to be stable? Try % an average of the two fixed points as a guess. % % The numer of total DNA Binding Sites u = zeros(INPUTS,1); DNATotal = 1; u(DNATotal, 1) = 1; % Indicies into the x, u, a1, a2 , b1 and b2 MCI = 1; CI = 2; CI2 = 3; CI2B = 4; MCRO = 5; CRO = 6; CRO2 = 7; CRO2B = 8; % Number of variables N = 8; % Number of inputs INPUTS = 1; % Initialize a1, a2 b1, b2, and u % We will fill in the nonzero values later a1 = zeros(N,N); a2 = zeros(N,N*N); b1 = zeros(N, INPUTS); b2 = zeros(N, INPUTS*N); u = zeros(INPUTS,1); I = eye(N); % Synth and Deg of MCI and MCRO [a1, a2] = produce(a1, a2, CI2B, MCI, 0.5); [a1, a2] = produce(a1, a2, CRO2B, MCRO, 0.5); [a1, a2] = deg(a1, a2, MCI, 1); [a1, a2] = deg(a1, a2, MCRO, 1); % Synth and deg of CI and CRO [a1, a2] = produce(a1, a2, MCI, CI, 1); [a1, a2] = produce(a1, a2, MCRO, CRO, 1); [a1, a2] = deg (a1, a2, CI, 1); [a1, a2] = deg (a1, a2, CRO, 1); % Dimerization of CI and CRO [a1, a2] = two2oneEq(a1, a2, CI, CI, CI2, 10, 1); [a1, a2] = two2oneEq(a1, a2, CRO, CRO, CRO2, 10, 1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Unfortunatly, the input needs to be handled specially % % In this case the binding to Cro2 and CI2 depends on the total ammount % of DNA binding sites. This is not controlled by the model so it is given % as an input. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% b2(CRO2B, (DNATotal-1)*INPUTS+CRO2) = 10; b2(CI2B, (DNATotal-1)*INPUTS+CI2) = 10; a2(CRO2B, (CRO2 -1)*N+CRO2B) = -10/2; a2(CRO2B, (CRO2B-1)*N+CRO2) = -10/2; a1(CRO2B, CRO2B) = - 1; b2(CRO2, (DNATotal-1)*INPUTS+CRO2) = -10; a2(CRO2, (CRO2 -1)*N+CRO2B) = 10/2; a2(CRO2, (CRO2B-1)*N+CRO2) = 10/2; a1(CRO2, CRO2B) = 1; a2(CRO2, (CRO2-1)*N+CI2B) = 10/2; a2(CRO2, (CI2B-1)*N+CRO2) = 10/2; a2(CRO2, (CRO2-1)*N+CI2B) = 10/2; a2(CRO2, (CI2B-1)*N+CRO2) = 10/2; a2(CRO2B, (CRO2-1)*N+CI2B) = -10/2; a2(CRO2B, (CI2B-1)*N + CRO2) = -10/2; a2(CI2B, (CI2-1)*N+CI2B ) = -10/2; a2(CI2B, (CI2B-1) *N+ CI2) = -10/2; a2(CI2B, (CI2-1) *N+CRO2B) = -10/2; a2(CI2B, (CRO2B-1)*N+CI2 ) = -10/2; a1(CI2, CI2B) = 1; a1(CI2B, CI2B) = - 1; b2(CI2, (DNATotal-1)*INPUTS+CI2) = -10; a2(CI2, (CI2-1) *N+CI2B) = 10/2; a2(CI2, (CI2B-1) *N+CI2) = 10/2; a2(CI2, (CI2-1) *N+CRO2B) = 10/2; a2(CI2, (CRO2B-1)*N+CI2) = 10/2; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % The Rate Law and the Jacobian of the Rate Law %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% dxdt = %% CHANGEME dxdx = %% CHANGEME %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Helper functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % k+ % A <----> B % k- function [a1, a2] = one2oneEq(a1, a2, from, to, fwdrate, revrate) a1(to, from) = fwdrate; a1(from, from) = -fwdrate; a1(from, to) = revrate; a1(to, to) = -revrate; % k+ % A ----> B % function [a1, a2] = one2one(a1, a2, from, to, fwdrate) a1(to, from) = fwdrate; a1(from, from) = -fwdrate; % k+ % A + B ----> C % function [a1, a2] = two2one(a1, a2, from1, from2, to, fwdrate) N = size(a2,1); if(from1 == from2) a2(to, N*(from1-1)+from1) = fwdrate; a2(from1, N*(from1-1)+from1) = -2 * fwdrate; else a2(to, N*(from1-1)+from2) = fwdrate/2; a2(to, N*(from2-1)+from1) = fwdrate/2; a2(from1, N*(from1-1)+from2) = -fwdrate/2; a2(from1, N*(from2-1)+from1) = -fwdrate/2; a2(from2, N*(from1-1)+from2) = -fwdrate/2; a2(from2, N*(from2-1)+from1) = -fwdrate/2; end % k+ % A + B <----> C % k- % function [a1, a2] = two2oneEq(a1, a2, from1, from2, to, fwdrate, revrate) N = size(a2,1); if(from1 == from2) a2(to, N*(from1-1)+from1) = fwdrate; a2(from1, N*(from1-1)+from1) = -2 * fwdrate; a1(from1, to) = 2 * revrate; a1(to, to) = - revrate; else a2(to, N*(from1-1)+from2) = fwdrate/2; a2(to, N*(from2-1)+from1) = fwdrate/2; a2(from1, N*(from1-1)+from2) = -fwdrate/2; a2(from1, N*(from2-1)+from1) = -fwdrate/2; a2(from2, N*(from1-1)+from2) = -fwdrate/2; a2(from2, N*(from2-1)+from1) = -fwdrate/2; a1(from1, to) = revrate; a1(from2, to) = revrate; a1(to, to) = - revrate; end % k+ % A ----> B + C % function [a1, a2] = one2two(a1, a2, from1, to1, to2, fwdrate) a1(to1, from1) = fwdrate; a1(to2, from1) = fwdrate; a1(from, from1) = -fwdrate; % k+ % A <----> B + C % k- function [a1, a2] = one2twoEq(a1, a2, from1, to1, to2, fwdrate) a1(to1, from1) = fwdrate; a1(to2, from1) = fwdrate; a1(from1, from1) = -fwdrate; a2(from1, N*(to1-1)+to2) = revrate/2; a2(from1, N*(to2-1)+to1) = revrate/2; a2(to1, N*(to1-1)+to2) = revrate/2; a2(to1, N*(to2-1)+to1) = revrate/2; a2(from1, N*(to1-1)+to2) = revrate/2; a2(from1, N*(to2-1)+to1) = revrate/2; a2(to2, N*(to1-1)+to2) = revrate/2; a2(to2, N*(to2-1)+to1) = revrate/2; % % First order decay % function [a1, a2] = deg(a1, a2, deg, rate) a1(deg, deg) = -rate; % % A produces B with first order kinetics % function [a1, a2] = produce(a1, a2, producer, product, rate) a1(product, producer) = rate;