% genereg.m - gene regulation example
% RMM, 4 Jul 07

% Parameters (see DynamicBehavior/geneticswitch.m)
genereg_params;				% get lacI parameters

alpha =  k_transcription_lacI;
alpha0 = k_transcription_leakage;
gamma = k_mRNA_degradation;
delta = k_protein_degradation;
beta = k_translation;
n = n_lacI;			
k = 1/(KM_lacI^n_lacI);

% Open loop transfer function
peq = (beta/delta) * ((alpha0)/gamma)
Gopen = zpk([], [-gamma, -delta], beta);

% Figure out the equilibrium point
peq = fzero(@(p) ((gamma*delta) / (alpha*beta) * p - alpha0/alpha) ...
  * (1 + k * p^n) - 1, peq)

% Closed loop transfer function (assumes Hill coefficient is n = 2)
sigma = 2*alpha*k*peq / (1 + k*peq^2)^2;
Gclosed = tf([beta], [1, gamma+delta, gamma*delta + sigma*beta]);

% Figure out a scale factor
scale = evalfr(Gopen, 0);

figure(1); clf; subplot(221);
W = logspace(-4, -1);
[Mopen, Popen] = bode(Gopen/scale, W);
[Mclosed, Pclosed] = bode(Gclosed/scale, W);
loglog(W, squeeze(Mopen), '--', W, squeeze(Mclosed), '-', ...
  'Linewidth', AM_bode_linewidth);
amaxis([1e-4 2e-2 0.01 2]);
xlabel('freq');
ylabel('mag');

lgh = legend('open loop', 'negative feedback', 'Location', 'south');
legend(lgh, 'boxoff');

amprint('genereg-bode.eps');
