function Freiburg2008_M1_Xten()
%% TCR Dimerization Model (steadystate/turnover/internalization/interfacepool/sparepool)
% written by Robert Gawlik, Freiburg iGEM team 2008
%%
clear all

%% Defining parameters
s = 0.1;                    %Turnover rate
Kon = 1;                    %Binding Rate
Koff = 0.2;                 %Dissociating Rate
Kdon = 1;                   %Dimerization Rate
Kdoff = 0.2;                %Dimer Dissociating Rate
Ka = 1;                     %Activation rate
Ki = 1;                     %Internalization Rate
N = 0.2;                    %Ligand input 
l = 1;                      %Ratio between areas of spare and interface pool
p = 1;                      %Exchange rate.
h = 2;                      %Kinetic order of the reactions.
 

%% Initial integration conditions
ti = 0:0.1:10;              %Time Vector for integrator

S0 = 1;                     %TCR in spare pool
T0 = 1;                     %TCR in interface pool
TA0 = 0;                    %Active TCR

%% Solving equations and plotting
[t,y]=ode23(@S_I_A,(ti),[S0 T0 TA0],[],N,h,Kon,Koff,Kdon,Kdoff,Ka,Ki,s,l,p);
figure;
plot(t,y,'linewidth',2);
title('TCR densities in time');
xlabel('time');
ylabel('TCR density');
legend('TCR in spare pool','TCR in interface pool','active TCR');

%% 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
    p=n;                         %using parameter p: increasing and solving ODEs in each loop cycle
    q=fix(n*10);
    disp(['solving for parameter =  ', num2str(n)]);
    [t,y] = ode23(@S_I_A,(ti),[S0 T0 TA0],[],N,h,Kon,Koff,Kdon,Kdoff,Ka,Ki,s,l,p);
    m(:,q)=y(:,3);               %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 = S_I_A(t,y,N,h,Kon,Koff,Kdon,Kdoff,Ka,Ki,s,l,p)
dydt = (zeros(size(y)));

S = y(1);                                 %TCR density in spare pool
T = y(2);                                 %TCR density in interface pool
TA = y(3);                                %Active TCR density
                                             
KG =2*Ka*(Kdon/(Kdoff+Ka));               %Combinating the rates
KI = Kon/Koff;

dydt(1) = -l*p*(S-T) + s*(1-S);                     %ODE for spare pool TCR
dydt(2) = p*(S-T) + s*(1-T) - KG*KI^h*T^h*N^h;      %ODE for interface pool TCR
dydt(3) = KG*KI^h*T^h*N^h  - Ki*TA;                 %ODE for active TCR
end