% genswitch_plot.m - run the Gardner/Collins genetic switch model
% RMM, 24 Jan 07

aminit;					% initialize MATLAB

% Run a simulation from an initial condition
sol1 = ode45(@genswitch, [0 10], [1 5]);

% Extract the equlibirum points
mu = 4; n = 2;			% switch parameters
Tfin = length(sol1.x);
eqpt = [sol1.y(1,Tfin), sol1.y(2,Tfin)];
eqpt(3) = fzero(@(x) mu/(1+x^2) - x, 2);

sol2 = ode45(@genswitch, [11 25], [sol1.y(1,Tfin)+2, sol1.y(2,Tfin)-2]);

% First plot out the curves that define the equilibria
u = 0:0.1:5;
f = mu ./ (1 + u.^n);

figure(1); clf; subplot(221);   hold on;
  amaxis([0 5 0 5]); box on;
  plot(u, f, '-', f, u, '--', 'LineWidth', AM_data_linewidth);
  lgh = legend('z1, f(z1)', 'z2, f(z2)'); legend(lgh, 'boxoff');
  plot([0, 3], [0, 3], 'k-', 'LineWidth', AM_ref_linewidth);
  plot(eqpt(1), eqpt(2), 'k.', eqpt(2), eqpt(1), 'k.', ...
    eqpt(3), eqpt(3), 'k.', 'MarkerSize', AM_data_markersize*3);
  xlabel('z1, f(z2)');
  ylabel('z2, f(z1)');

amprint('genswitch-eqpts.eps');

% Time traces
figure(3); clf; subplot(221);
  pl = plot(sol1.x, sol1.y(1,:), '-', sol1.x, sol1.y(2,:), '--');
  set(pl, 'LineWidth', AM_data_linewidth);
  hold on;
  pl = plot([sol1.x(Tfin), sol1.x(Tfin)+1], ...
       [sol1.y(1,Tfin), sol2.y(1,1)], 'o:', ...
       [sol1.x(Tfin), sol1.x(Tfin)+1], [sol1.y(2,Tfin), sol2.y(2,1)], 'o:');
  set(pl, 'LineWidth', AM_data_linewidth, 'MarkerSize', AM_data_markersize);
  pl = plot(sol2.x, sol2.y(1,:), '-', sol2.x, sol2.y(2,:), '--');
  set(pl, 'LineWidth', AM_data_linewidth);
  amaxis([0 25 0 5]);

  xlabel('Time {\itt} [scaled]');
  ylabel('Protein concentrations [scaled]');
  legh = legend('z1 (A)', 'z2 (B)', 'Orientation', 'horizontal');
  legend(legh, 'boxoff');

amprint('genswitch-time.eps');

% Phase portrait
figure(2); clf; subplot(221);
  amaxis([0 5 0 5]); set(gca, 'DataAspectRatio', [1 1 1]);
  amphaseplot('genswitch', ...
    'timepts', [0.2 0.6 1.2], nan, ...
    boxgrid([0, 5, 6], [0, 5, 6]), 10 ...
  );

  % Add the stable equilibrium points
  hold on;
  plot(eqpt(1), eqpt(2), 'k.', eqpt(2), eqpt(1), 'k.', ...
    eqpt(3), eqpt(3), 'k.', 'MarkerSize', AM_data_markersize*3);

  xlabel('Protein A [scaled]');
  ylabel('Protein B [scaled]', 'Rotation', 90);

amprint('genswitch-phase.eps');

