%% 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 the Huang & Ferrell model of the MAPK cascade.
% model = loadModel('ThreeStepMAPK');
model = loadModel('../../models/EGFR model/Apgar2008');

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

%% 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) + 4, 100);
[yss, xss, val] = doseResponse3(model, u, x(tF)); 

%% Plot Results
% Again, we plot the dose response curves for the system outputs.  In this
% case, the outputs are the active MAPKKK, MAPKK, and MAPK species.
figure;
semilogx(u, yss,  '-');
legend(model.yNames);
xlabel('Ligand Concentration (L)');
ylabel('Ouput Concentration');
title('Dose Response Curve for Three-Step Model (Outputs Only)');
