function Freiburg2008_MSRS()
%% Modular Synthetic Receptor System model
% written by Robert Gawlik, Freiburg iGEM team 2008

%% 
clear all

%% Defining parameters
s1 = 0.1;          %turnover rate for R1
s2 = 0.1;          %turnover rate for R2
K1 = 3;            %binding rate 1
k1 = 0.1;          %dissociating rate 1
K2 = 3;            %binding rate 2
k2 = 0.1;          %dissociating rate 2
K3 = 2;             %R1R2 dimerization rate 1
k3 = 0.1;           %R1R2 dimer dissociating rate 1
K4 = 2;             %R1R2 dimerization rate 2
k4 = 0.1;           %R1R2 dimer dissociating rate 2
K5 = 2;             %R1R1 dimerization rate 
k5 = 0.1;           %R1R1 dimer dissociating rate 
K6 = 2;             %R2R2 dimerization rate 
k6 = 0.1;           %R2R2 dimer dissociating rate 

Ka = 1;            %activation rate
Kda = 0.05 ;       %deactivation rate

%% initial integration conditions
ti=0:0.1:10;

R10 = 1;            %receptor1
R20 = 1;            %receptor2
NN0 = 1;            %ligand
R1NN0 = 0;          %monomer 1
R2NN0 = 0;          %monomer 2
R1NNR20 = 0;        %dimer 1
R1NNR10 = 0;        %dimer 2
R2NNR20 = 0;        %dimer 3
A0 = 0;             %active protein
%% solving and plotting
[t,y]=ode45(@MSRS_ODEs,(ti),[R10 R20 NN0 R1NN0 R2NN0 R1NNR20 R1NNR10 R2NNR20 A0],[],s1,s2,K1,k1,K2,k2,K3,k3,K4,k4,K5,k5,K6,k6,Ka,Kda);
figure;
plot(t,y(:,1),'blue');
hold on
plot(t,y(:,2),'-.blue');
hold on
plot(t,y(:,3),'green');
hold on
plot(t,y(:,4),'black');
hold on
plot(t,y(:,5),'-.black');
hold on
plot(t,y(:,6),'cyan','linewidth',2);
hold on
plot(t,y(:,7),'yellow');
hold on
plot(t,y(:,8),'-.yellow');
hold on
plot(t,y(:,9),'red','linewidth',2);
hold on
legend('free R1','free R2','ligand','R1-ligand monomer','R2-ligand monomer ','receptor dimer R1R2','receptor dimer R1R1','receptor dimer R2R2','active protein');
xlabel('time');
ylabel('densities of involved quantities');
title('Modular synthetic receptor system activity')

%% 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:3                  %parameter range 0.1 to 3
    NN0=n;                       %using parameter NN0: increasing and solving ODEs in each loop cycle
    q=fix(n*10);
    disp(['solving for parameter =  ', num2str(n)]);
    [t,y]=ode45(@MSRS_ODEs,(ti),[R10 R20 NN0 R1NN0 R2NN0 R1NNR20 R1NNR10 R2NNR20 A0],[],s1,s2,K1,k1,K2,k2,K3,k3,K4,k4,K5,k5,K6,k6,Ka,Kda);
    m(:,q)=y(:,9);               %taking each cycle ODE solution of quantity of interest for plotting
    subplot(5,6,q);
    plot(t,y);
end   
figure
surf(m','meshstyle','row');
xlabel('time');
ylabel('parameter');
zlabel('density');
end
%% ODEs
function dydt = MSRS_ODEs(t,y,s1,s2,K1,k1,K2,k2,K3,k3,K4,k4,K5,k5,K6,k6,Ka,Kda)
    dydt=zeros(size(y));

    R1 = y(1);           %Receptor 1
    R2 = y(2);           %Receptor 2
    NN = y(3);           %Ligand
    R1NN = y(4);         %Receptor-Ligand complex 1
    R2NN = y(5);         %Receptor-Ligand complex 2
    R1NNR2 = y(6);       %Receptor-Ligand dimer
    R1NNR1 = y(7);       %Receptor-Ligand dimer Err /no act
    R2NNR2 = y(8);       %Receptor-Ligand dimer Err /no act
    A = y(9);            %active receptor

    dydt(1) = s1*(1-R1) + k1*R1NN - K1*R1*NN + k4*R1NNR2 - K4*R1*R2NN + k5*R1NNR1 - K5*R1*R1NN +Kda*A^2;
    dydt(2) = s2*(1-R2) + k2*R2NN - K2*R2*NN + k3*R1NNR2 - K3*R2*R1NN + k6*R2NNR2 - K6*R2*R2NN +Kda*A^2;
    dydt(3) = -K1*R1*NN + k1*R1NN - K2*R2*NN + k2*R2NN;
    dydt(4) = K1*R1*NN - k1*R1NN + k3*R1NNR2 - K3*R2*R1NN + k5*R1NNR1 - K5*R1*R1NN + Kda*A^2;
    dydt(5) = K2*R2*NN - k2*R2NN + k4*R1NNR2 - K4*R1*R2NN + k6*R2NNR2 - K6*R2*R2NN + Kda*A^2;
    dydt(6) = K3*R2*R1NN - k3*R1NNR2 + K4*R1*R2NN - k4*R1NNR2 - Ka*R1NNR2;
    dydt(7) = K5*R1*R1NN - k5*R1NNR1;
    dydt(8) = K6*R2*R2NN - k6*R2NNR2;
    dydt(9) = Ka*R1NNR2 - 2*Kda*A^2;
end