function Freiburg2008_M1_Basic()
%% Model 1
% written by Robert Gawlik, Freiburg iGEM team 2008
%%
clear all
%% Defining parameters:
s = 0.1;                 %turnover rate
kon = 1;                %TCR NIP binding rate
koff = 0.2;               %TCR NIP dissociating rate
kdon = 1;               %TCR-NIP dimerization rate
kdoff = 0.2;              %TCR-NIP dissociating rate
ka = 1;                 %activation rate
ki = 0.8;                %internalisation rate
N = 0.2;                 %NIP amount
h = 2;                  %kinetic order 

%% initial integration conditions
ti = 0:0.1:10;
T0 = 1;
TA0 = 0;

%% integrating and plotting
[t,y]=ode45(@T_A,(ti),[T0 TA0],[],N,h,s,kon,koff,kdon,kdoff,ka,ki);
figure;
plot(t,y(:,1));
hold on;
plot(t,y(:,2),'red','linewidth',2);
xlabel('time')
ylabel('TCR densities');
title('Active TCR in time');
legend('free TCR','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
    s=n;                         %using parameter s: increasing and solving ODEs in each loop cycle
    q=fix(n*10);
    disp(['solving for parameter =  ', num2str(n)]);
    [t,y] = ode45(@T_A,(ti),[T0 TA0],[],N,h,s,kon,koff,kdon,kdoff,ka,ki);
    m(:,q)=y(:,2);               %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_A(t,y,N,h,s,kon,koff,kdon,kdoff,ka,ki)
dydt=zeros(size(y));

T = y(1);                                                        %membrane TCR 
TA= y(2);                                                        %active TCR

KG = 2*ka*(kdon/(kdoff+ka)); 
KI = kon/koff;

dydt(1) = s*(1-T) - KG*KI^h*T^h*N^h;                             %membrane TCR building
dydt(2) = KG*KI^h*T^h*N^h - ki*TA;                               %activation of TCR
end