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