% INDUCIBLE MODEL LuxR-LuxI RESOLUTION:

% Parameters
time_effect =  0.017;      % doubling time effect [1/min]
NA = 6.023e23;                          % Avogadro's number: # particles/mol
Vcell = 1.1e-9;                   % typical volume of E. coli (microliters/cell). Source: Bionumbers
Vext = 1e-3;                     % culture medium volume [ul] in a microfluidics device
Ccells_OD600_1=8e5;               % cells/microliter for OD=1. 
                                        % Source: Agilent, E. coli Cell Culture Concentration 
                                        % from OD600 Calculator
Ncell_0 = 0.1*Ccells_OD600_1*Vext;                                  
factor_units=1e15/(Vcell*NA); % Conversion factor from number of particles to
                                  % concentration (1 nMolar=nanomols/liter)
                                  % 1 particle/cell x 1/Vcell cell/microliter 
                                  % x 1/nA mol/particle x
                                  % 1e6 microliter/liter x 1e9 nanomol/l

pN_luxI = 17;                 % plasmid number of luxI in pBR322 (15-20 copies/cell)
pN_luxR = 10;                 % plasmid number of luxI in pACYC184 (10 copies/cell)

% LuxI
alphaI = 0.17;              % leakage of repressor
dmI = time_effect + 0.23;                % degradation rate of LuxI mRNA [1/min]. Typical mRNA degradation rate for E coli
dI = time_effect + 0.01;                  % degradation rate of LuxI [1/min]
pI = 3.09;                     % translation rate of LuxI #mRNA [1/min]. b*dmI from Weber.  [3.0928 - 6.1856]
kI = 1.03;                      % average transcription rate of LuxI [1/min] from calculator [1.0309 - 10.3093] from our rates calculator based on Alberts, and
                                % Ref 2: 50/b from Weber
% LuxR
alphaR = 0.01;             % ratio between unactivated vs activated of expresion LuxR
dmR = time_effect + 0.23;                % degradation rate of LuxR mRNA [1/min]
dR = time_effect + log(2)/5;                 % degradation rate of LuxR [1/min]
pR =  2.38; % b*dmR;                  % translation rate of LuxR [1/min] [2.381 - 4.7619] from our rates calculator, Alberts
kR = 0.79; % 200/b;                 % transcription rate of LuxR [1/min] [0.79365 - 7.9365] from our rates calculator
cR = kR*pN_luxR;  % constitutive expression of LuxR (always ConstR = initCond of plasmids in complete model)
% Monomer
kd1 = 10; % 100;                  % dissociation constant of R to A [nM], Urbanowski etal. 2004
k_1 = 1;% 10;                    % unbinding rate LuxR to AHL [1/min]

dRA = time_effect + log(2)/5;   % degradation rate of (LuxR.A) [1/min]. Buchler et al. 2004 Monomer half-life is just few minutes.

% Dimer
kd2 = 10;%20;                   % or 10 dissociation cte (LuxR.A) [nM], Buchler et al. 2003
                            %Koren, R. & Hammes, G. G. (1976) Biochemistry 15, 1165–1171.
                            %Northrup, S. H. & Erickson, H. P. (1992) Proc. Natl. Acad. Sci. USA 89, 3338–3342.
k_2 = 0.1;                     % dissociation rate dimer (LuxR.A)2 [1/min]


kdLux = 10;           % dissociation cte (LuxR.A)2 to promoter [nM], Bucler et al [1 1000]nM
k_Lux = 10;                  % dissociation rate(LuxR.A)2 to Lux promoter [1/min]
kLux = k_Lux/kdLux;               % binding rate LuxR to AHL [1/min]
dRA2 = time_effect;         % Ron Weiss et al. A synthetic multicellular system for programmed pattern formation
                            % Buchler 2004 Degradation rate for (LuxR.AHL)_2. Corresponding to dilution if cell half-life = tau minutes 

D = 2;                    % kinetic rate of AHL external transport [1/min] across the cell membrane, calculated
dA = time_effect + 0.0002;                 %[0.05 0.03 min^-1]Degradation from Bionumbers online      
                            %0.0164, Degradation rate for AHL. From Kaufmann etal. 2005. Similar to You etal. Nature, 2004 
dAe = 0.0002;  
Vc = Vcell/Vext;
N = 240; % cells

param = [pI, kI, pN_luxI, dmI, kdLux, alphaI, dI, k_1, kd1, pR, cR, dmR, k_2, kd2, dRA, dR, dRA2, D, Vc, dA, dAe, N];

%% SIMULATION

AHL = 100; % en nM
AHL0 = AHL*NA*1e-9*1e-9  ; % AHL nM x 1e-9M/1nM x N Avogadro x Vol (1e-9L)
% Initial Conditions
InitC = [929, 505, 0, 0, AHL0];

time = 1:0.1:300;
options = odeset('AbsTol',1e-6,'RelTol',1e-3);
[t,n] = ode15s(@(t,x) mi_lux(x,param),time,InitC,options);

% Graphs
fig = figure('Color','w');
hold on
subplot(2,3,1)
plot(t,n(:,1),'r')
xlabel('Time (min)','FontName','Lato')
ylabel('PoI (molec)','FontName','Lato')
title('PoI vs Time','FontName','Lato','FontSize',14)
hold on
subplot(2,3,2)
plot(t,n(:,2),'r')
xlabel('Time (min)','FontName','Lato')
ylabel('LuxR (molec)','FontName','Lato')
title('LuxR vs Time','FontName','Lato','FontSize',14)
hold on
subplot(2,3,3)
plot(t,n(:,3),'k')
xlabel('Time (min)','FontName','Lato')
ylabel('(LuxR.AHL)_2 (molec)','FontName','Lato')
title('(LuxR.AHL)_2 vs Time','FontName','Lato','FontSize',14)
hold on
subplot(2,3,4)
plot(t,n(:,4),'r')
xlabel('Time (min)','FontName','Lato')
ylabel('AHL_i (molec)','FontName','Lato')
title('AHL_i vs Time','FontName','Lato','FontSize',14)
hold on
subplot(2,3,5)
plot(t,n(:,5),'b')
xlabel('Time (min)','FontName','Lato')
ylabel('AHL_e (molec)','FontName','Lato')
title('AHL_e  vs Time','FontName','Lato','FontSize',14)
hold on

%% ANALYSIS
load data_exp_ahl.mat

AHL = logspace(-2,3,100); % nM
AHL0 = AHL*NA*1e-9*1e-9  ; % AHL nM x 1e-9M/1nM x N Avogadro x Vol (1e-9L)
param = [pI, kI, pN_luxI, dmI, kdLux, alphaI, dI, k_1, kd1, pR, cR, dmR,...
        k_2, kd2, dRA, dR, dRA2, D, Vc, dA, dAe, N];
PoI = [];  
for i=1:length(AHL)
    % Initial Conditions
    InitC = [929, 505, 0, 0, AHL0(i)];

    time = 1:0.1:300;
    options = odeset('AbsTol',1e-6,'RelTol',1e-3);
    [t,n] = ode15s(@(t,x) mi_lux(x,param),time,InitC,options);
    PoI = [PoI n(end,1)];
end

% Graphs
blue = [68 146 198]/255;
red = [246 80 88]/255;
yellow = [255 198 0]/255;

fig1 = figure('Color','w');
PoI = PoI/max(PoI);
semilogx(AHL,PoI,'Color',blue,'LineWidth',1.8)
hold on
semilogx(data_exp(:,1),data_exp(:,2),'MarkerSize',6,...
    'MarkerFaceColor',red, 'MarkerEdgeColor',[0 0 0],...
    'Marker','o', 'LineStyle','none')
xlabel('AHL concentration (nm)','FontName','Lato')
ylabel('GFP normalised fluorescence (a.u)','FontName','Lato')
title('GFP  vs AHL','FontName','Lato','FontSize',16)
legend({'Fitted model','Experimental data'},'Location','northeast')
grid on