function Freiburg2008_M_2()
%% Dimerization model 2 
% written by Robert Gawlik, Freiburg iGEM team 2008
%%
clear all
%% Defining parameters
s = 0.1;
Kon = 3;                %receptor ligand binding rate
Koff = 0.1;             %receptor ligand dissociating rate
Kdon = 1;               %dimerization rate
Kdoff = 0.2;            %dimer dissociating rate
Ka = 1;                 %activation rate
Ki = 0.8;               %internalisation rate 
%% initial integrating contitions
ti=0:0.1:10;

R0 = 1;         %receptor
NN0 = 0.5;      %bivalent ligand
RNN0 = 0;       %monomer
RNNR0 = 0;      %dimer
A0 = 0;         %activ
%% solving and plotting
[t,y]=ode23(@R_A,(ti),[R0 NN0 RNN0 RNNR0 A0],[],s,Kon,Koff,Kdon,Kdoff,Ka,Ki);
figure;
plot(t,y);
legend('Receptor','NIP','Receptor-NIP','Receptor dimer','active receptor');

figure
plot(t,y(:,1),'blue');
hold on
plot(t,y(:,3),'black');
hold on
plot(t,y(:,4),'green');
hold on
plot(t,y(:,5),'red','linewidth',2);
legend('free receptor','receptor-NIP monomer','receptor dimer','active receptor');
%legend('free receptor','receptor dimer','active receptor');
xlabel('time');
ylabel('densities');
title('Receptor densities in time');

%% loop for parameter analysis
%within this cell a parameter of interest can be set to the loop variable n
%to obtain several solutions with different values for the parameter.

figure
for n=0.1:0.1:2                  %parameter range 0.1 to 2
    Ki=n;                        %using parameter Ki: increasing and solving ODEs in each loop cycle
    q=fix(n*10);
    disp(['solving for parameter =  ', num2str(n)]);
    [t,y]=ode23(@R_A,(ti),[R0 NN0 RNN0 RNNR0 A0],[],s,Kon,Koff,Kdon,Kdoff,Ka,Ki);
    m(:,q)=y(:,5);               %taking each cycle ODE solution of quantity of interest for plotting
    subplot(5,4,q);
    plot(t,y);
end   
figure
surf(m','meshstyle','row');
xlabel('time');
ylabel('parameter');
zlabel('density');
end

%% ODEs
function dydt = R_A(t,y,s,Kon,Koff,Kdon,Kdoff,Ka,Ki)
    dydt=zeros(size(y));

    R = y(1);           %Receptor
    NN = y(2);          %Ligand
    RNN = y(3);         %Receptor-Ligand complex
    RNNR = y(4);        %Receptor-Ligand dimer
    A = y(5);           %active receptor

    dydt(1) = s*(1-R) + Koff*RNN - Kon*R*NN + Kdoff*RNNR - Kdon*RNN*R; 
    dydt(2) = Koff*RNN - Kon*R*NN;
    dydt(3) = Kon*R*NN + Kdoff*RNNR - Koff*RNN - Kdon*R*RNN; 
    dydt(4) = Kdon*RNN*R - Kdoff*RNNR - Ka*RNNR;
    dydt(5) = 2*Ka*RNNR - Ki*A;
end