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