%% Demo Dose Response Command
%
% <<../../../images/LogoSmall.jpg>>
%
% This scripts demonstrates how to use the |doseResponse| command.  This 
% command is used to find the steady state concentrations for a model
% given a particular input u.  Mathematically this means solving for
% x that is consistent with the initial conditions such that:
%
% $$\dot{x} = A_1 x + A_2 x \otimes  x + B_1 u + B_2 x \otimes u + K = 0$$
%
% We solve this equation using Newton's method which allows the method to
% find both stable and unstable fixed points.
%

%% Load the model
% First we load the model using the |loadModel| command.  In this case it
% is a model of a distributive MAPK phosphorylation dephosphophorylation
% loop.

tF = 1e3;
u0 = 1e-2;
% model = loadModel('OneStepMAPK');
model = loadModel('ThreeStepMAPK');

% tF = 1e5;
% u0 = 1e-15;
% model = loadModel('Apgar2008');

% tF = 1e5;
% u0 = 1e-4;
% model = loadModel('ferrell');

exp       = constructExperiment(tF, u0);
[t, y, x] = simulate(model, exp, struct('rel_tol', 1e-9, 'abs_tol', 1e-12));

plot(t, y(t)), pause



% u   = logspace(log10(u0), log10(u0)+5, 1e2);
% yss = zeros(length(u), model.nY);
% xss = zeros(length(u), model.nX);
% for i = 1:length(u)
%     i
%     exp = constructExperiment(tF, u(i));
%     [ti, yi, xi] = simulate(model, exp, struct('rel_tol', 1e-9, 'abs_tol', 1e-12));
%     yss(i,:) = yi(t(end));
%     xss(i,:) = xi(t(end));
%     figure(1)
%     semilogx(u(1:i), yss(1:i,:),  '--');
%     set(gca, 'XLim', [u(1) u(end)]);
%     drawnow
% end

%% Compute the dose response curve
% We solve for the steady state concentration of the system over a range
% of kinase concentrations u.  The output is the steady state
% concentration as well as the residual of the steady state equation.
%
u               = logspace(log10(u0), log10(u0)+5, 1e2);
% [yss_conv, xss_conv, val] = doseResponse3(model, u, x(tF)); 
[uss_conv, yss_conv, xss_conv, val] = doseResponse(model, u0, u0*1e4, x(tF)); 

%% Plot Results
% We plot the dose response curves for the System outputs, in this case
% the unphosphorylated substrate (S), the singly phosphorylated substrate
% (Sp) and the doubly phosphorylated substrate (Spp).  The second plot
% shows the dose response curves of all of the model species. 
%
figure(1)
semilogx(uss_conv, yss_conv,  '.-');
legend(model.yNames);
xlabel('Ligand Concentration (L)');
ylabel('Ouput Concentration');
title('Dose Response Curve for One-Level MAPK Model (Outputs Only)');
