function Freiburg2008_M1_Corrected()
%% Dimerization model 1 with correction of Bachmann error
% written by Robert Gawlik, Freiburg iGEM team 2008
%%
clear all

%% Defining constants

s = 0.1;                  %turnover rate
kon = 3;                  %TCR-NIP binding rate
koff = 0.1;               %TCR-NIP dissociating rate
kdon = 1;                 %dimerizing rate
kdoff = 0.2;              %dimer dissociating rate
ka = 1;                   %TCR activation rate
ki = 0.8;                 %internalization rate

%% initial integration conditions

ti = 0:0.1:10;     %time vector for integration

N0 = 1;            %ligand 

T0 = 1;           %TCR in membrane
TN0 = 0;          %TCR-NIP complex
TND0 = 0;         %Dimer
TA0 = 0;          %active TCR

%% integration and plotting

[t,y] = ode45(@T_TT_A,(ti),[N0 T0 TN0 TND0 TA0],[],s,kon,koff,kdon,kdoff,ka,ki);
figure
plot(t,y);
legend('NIP','free TCR','TCR-NIP complex','TCR-dimer','active TCR');

figure
plot(t,y(:,1),'green');
hold on
plot(t,y(:,2),'blue');
hold on
plot(t,y(:,4),'black');
hold on
plot(t,y(:,5),'red','linewidth',2);
legend('free TCR','NIP','dimer','active TCR');
xlabel('time');
ylabel('densities');
title('Free TCR, NIP and active TCR 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
    N0=n;                        %using parameter N0: increasing and solving ODEs in each loop cycle
    q=fix(n*10);
    disp(['solving for parameter =  ', num2str(n)]);
    [t,y] = ode23(@T_TT_A,(ti),[N0 T0 TN0 TND0 TA0],[],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 = T_TT_A(t,y,s,kon,koff,kdon,kdoff,ka,ki)
dydt = (zeros(size(y)));


N = y(1);           %NIP
T = y(2);           %TCR in the membrane
TN = y(3);          %TCR-NIP complex
TND = y(4);         %TCR-dimer 
TA = y(5);          %active TCR


dydt(1) = koff*TN - kon*T*N;                                 %ligand 
dydt(2) = s*(1-T) + koff*TN - kon*T*N;                       %receptor
dydt(3) = kon*T*N - koff*TN - 2*kdon*TN.^2 + 2*kdoff*TND;    %monomer
dydt(4) = kdon*TN.^2 - kdoff*TND - ka*TND;                   %dimer
dydt(5) = 2*ka*TND - ki*TA;                                  %active ones
end