%iGEM 2008 :MB

close all; 
clear all;
clc;
warning off MATLAB:divideByZero

global K_1 K_2 K_3 K_4 K_5 K_6 K_7 K_8 delta1 delta2;



disp(' ');
disp('__________________________________________________________________________');
disp('iGEM 2008 : Inducible Promoter Model 1');
disp('__________________________________________________________________________');
disp(' ');

tmax = input('t_max of the simulation? (Default is 10) ');
if isempty(tmax)
    tmax = 10;
end

K_1 = input('K_1 ? (Default is 1) ');
if isempty(K_1)
    K_1 = 1;
end
K_2 = input('K_2 ? (Default is 1) ');
if isempty(K_2)
    K_2 = 1;
end
K_3 = input('K_3 ? (Default is 1) ');
if isempty(K_3)
    K_3 = 1;
end
K_4 = input('K_4 ? (Default is 1) ');
if isempty(K_4)
    K_4 = 1;
end
K_5 = input('K_5 ? (Default is 1) ');
if isempty(K_5)
    K_5 = 1;
end
K_6 = input('K_6 ? (Default is 1) ');
if isempty(K_6)
    K_6 = 1;
end
K_7 = input('K_7 ? (Default is 1) ');
if isempty(K_7)
    K_7 = 1;
end
K_8 = input('K_8 ? (Default is 1) ');
if isempty(K_8)
    K_8 = 1;
end

delta1 = input('Degradation rate of LacI ? (Default is 1) ');
if isempty(delta1)
    delta1 = 1;
end

delta2 = input('Degradation rate of GFP ? (Default is 1) ');
if isempty(delta2)
    delta2 = 1;
end



disp(' ');
disp('__________________________________________________________________________');
disp(' Initial Conditions');
disp('__________________________________________________________________________');
disp(' ');


x0(1)=K_1/delta1;
alpha = input('Initial Concentrationof IPTG? (Default is 1) ');
if isempty(alpha)
    alpha= 1;
end
x0(2)=alpha;
x0(3)=0;

alpha = input('Initial Concentrationof of Promoter? (Default is 1) ');
if isempty(alpha)
    alpha= 1;
end
P_0=alpha;

K_alpha=K_2/K_3;
K_beta=K_4/K_5;
K_gamma=K_6/K_7;

x0(4)= P_0/(1+K_beta*x0(1)*x0(1));
x0(5)= x0(4)*K_beta*x0(1)*x0(1);
x0(6)=K_8*x0(4)/delta2;

% x(1); % LacI
% x(2); % IPTG
% x(3); % LacI2-IPTG
% x(4); % P
% x(5); % P-LacI2
% x(6); % GFP

disp(' ');
disp('__________________________________________________________________________');
disp(' Preliminary Predictions');
str=[' K_alpha= ' num2str(K_alpha) '   K_beta= ' num2str(K_beta) '   K_gamma= ' num2str(K_gamma)];
disp(str);

a=x0(1)*x0(1);
b=(a*K_2+K_3)/P_0+K_6;
X_1=1/(a*K_beta+1);
X_infinity=1/(1+a*K_alpha/K_gamma);

if (K_gamma==1)
    str=[' Root1 of Q = ' num2str(X_1) '; No Root2 of Q ;  Root-infinity =' num2str(X_infinity)];
else
    X_2=b/(K_6-K_7);
    str=[' Root1 of Q = ' num2str(X_1) '; Root2 of Q =' num2str(X_2)   '; Root-infinity =' num2str(X_infinity)];
end
disp(str);

if (K_beta>K_alpha/K_gamma)
    str2=['The limit GFP concentration is higher than the initial concentration '];
else
    str2=['The limit GFP concentration is lower than the initial concentration '];
end
disp(str2);
disp('__________________________________________________________________________');
disp(' ');

%__________________________________________________________________________
%__________________________________________________________________________
% Display Options: User chooses what plots they want to see

disp(' ');
disp('__________________________________________________________________________');
disp(' ');
disp('You can now choose what plots to display');
display_option = input('Do You want to see all the plots or customize the display ? (0=All else=Custom)');
if isempty(display_option)
    display_option = 0;
end

display_option_1=1;
display_option_2=1;
display_option_3=1;
display_option_4=1;
display_option_5=1;
display_option_6=1;


if (display_option>0)
    alpha = input('Do You want to see the Evolution of the concentration of LacI? (0=N else Y)');
    if isempty(alpha)
        alpha = 1;
    end
    if (alpha>0)
        alpha=1;
    end
    display_option_1=alpha;

    alpha = input('Do You want to see the Evolution of the concentration of IPTG? (0=N else Y)');
    if isempty(alpha)
        alpha = 1;
    end
    if (alpha>0)
        alpha=1;
    end
    display_option_2=alpha;
    
    alpha = input('Do You want to see the Evolution of the concentration of LacI2-IPTG Complex? (0=N else Y)');
    if isempty(alpha)
        alpha = 1;
    end
    if (alpha>0)
        alpha=1;
    end
    display_option_3=alpha;
%     
    alpha = input('Do You want to see the Evolution of the concentration of free Promoters else Y)');
    if isempty(alpha)
        alpha = 1;
    end
    if (alpha>0)
        alpha=1;
    end
    display_option_4=alpha;


    alpha = input('Do You want to see the Evolution of the concentration of the P-LacI2 Complex? (0=N else Y)');
    if isempty(alpha)
        alpha = 1;
    end
    if (alpha>0)
        alpha=1;
    end
    display_option_5=alpha;
    
    
    alpha = input('Do You want to see the Evolution of the concentration of FP? (0=N else Y)');
    if isempty(alpha)
        alpha = 1;
    end
    if (alpha>0)
        alpha=1;
    end
    display_option_6=alpha;
end

[t,x] = ode45('InduciblePromoter_ODE',[0 tmax],x0);

%__________________________________________________________________________
%__________________________________________________________________________
% Display of the Interesting Quantities

if (display_option_1>0)
    figure;
    maxi=0;
    plot(t, x(:,1));
    maxi=max(x(:,1));
    if (maxi==0)
        maxi=0.1;
    end
    
    axis([0 tmax 0 maxi]);
    xlabel('time');
    ylabel('Concentration of the concentration of LacI');
    title(' Evolution of the Concentration of the concentration of LacI ');
    hold off;
end

if (display_option_2>0)
    figure;
    maxi=0;
    plot(t, x(:,2));
    maxi=max(x(:,2));
    if (maxi==0)
        maxi=0.1;
    end
    

    axis([0 tmax 0 maxi]);
    xlabel('time');
    ylabel('Concentration of the concentration of IPTG');
    title(' Evolution of the Concentration of the concentration of IPTG');
    hold off;
end

if (display_option_3>0)
    figure;
    maxi=0;
    plot(t, x(:,3));
    maxi=max(x(:,3));
    if (maxi==0)
        maxi=0.1;
    end
    
    axis([0 tmax 0 maxi]);
    xlabel('time');
    ylabel('Concentration of the concentration of the LacI2-IPTG Complex');
    title(' Evolution of the Concentration of the concentration of the LacI-IPTG Complex');
    hold off;
end

if (display_option_4>0)
    figure;
    maxi=0;
    plot(t, x(:,4));
    maxi=max(x(:,4));
    if (maxi==0)
        maxi=0.1;
    end
    
    axis([0 tmax 0 maxi]);
    xlabel('time');
    ylabel('Concentration of the concentration of free Promoters');
    title(' Evolution of the Concentration of the concentration of free Promoters');
    hold off;
end

if (display_option_5>0)
    figure;
    maxi=0;
    plot(t, x(:,5));
    maxi=max(x(:,5));
    if (maxi==0)
        maxi=0.1;
    end
    
    axis([0 tmax 0 maxi]);
    xlabel('time');
    ylabel('Concentration of the concentration of the P-LacI Complex');
    title(' Evolution of the Concentration of the concentration of the P-LacI Complex');
    hold off;
end

if (display_option_6>0)
    figure;
    maxi=0;
    plot(t, x(:,6));
    maxi=max(x(:,6));
    if (maxi==0)
        maxi=0.1;
    end
    
    axis([0 tmax 0 maxi]);
    xlabel('time');
    ylabel('Concentration of Fluorescent Protein FP');
    title(' Evolution of the Concentration of Fluorescent Protein FP');
    hold off;
end

clear all
%__________________________________________________________________________
%__________________________________________________________________________