% Delay differential equation for the lac operon
% RMM, 20 Jan 07
%
% Based on Yildirim et al, Chaos 2004

function xder = lac(t, x, xlag)

xdel_M = xlag(:,1);
xdel_B = xlag(:,2);

% Parameters (from YSHM04; converted to uM)
mu_bar = 3.03e-2;			% dilution rate
alpha_M = 0.997;   gamma_M = 0.411;	% mRNA growth/decay rates
n = 2;					% cooperativity
K = 7200;				% basal transcription constant
K_1 = 2.52e-2;				% Hill constant for lacI binding
alpha_B = 1.66e-2; gamma_B = 8.33e-4;	% beta-gal growth/decay
alpha_A = 1.76e4;  gamma_A = 1.35e-2;	% allolactose growth/decay
beta_A = 2.15e4;
K_L = 970;	   K_A = 1950;		% allolactose Michaelis Menten terms
L = 100;				% external lactose concentration

tau_M = 0.1; tau_B = 2;			% time delays

% mRNA concentration
dM = alpha_M * ...
       (1 + K_1 * (exp(mu_bar*tau_M) * xdel_M(3))^n ) / ...
       (K + K_1 * (exp(mu_bar*tau_M) * xdel_M(3))^n ) - gamma_M * x(1);

% B-gal concentration
dB = alpha_B * exp(-mu_bar*tau_B) * xdel_B(1) - gamma_B * x(2);

% Allolactose concentration
dA = alpha_A * x(2) * L / (K_L + L) - ...
     beta_A * x(2) * x(3) / (K_A + x(3)) - gamma_A * x(3);

xder = [dM; dB; dA];