%% 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 concenetrations for a model
% given a particular input u.  Mathematically this means solving for
% x that is consitent 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.
%

%% Add the KroneckerBio toolbox to the path 

addpath ../../src/

%% Load the model
% First we load the model using the |loadModel| command.  In this case it
% is a model of a a distributive MAPK phosphorylation dephosphophorylation
% loop.
%model = loadModel('OneStep');
model = loadModel('..\..\models\EGFR model\Apgar2008');                
%% 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 residulat of the steady state equation.
%
u               = logspace(-14, -7, 100);
[yss, xss, val] = doseResponse(model, u, model.ic); 

%% 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 phosphorlated substrate (Spp).  The second plot
% shows the dose response curves of all of the model species. 
%
figure;
semilogx(u, yss*diag(1./max(yss)),  '-');
legend(model.yNames);
xlabel('Ligand Concentration (L)');
ylabel('Ouput Concentration');
title('Dose Response Curve for One-Step Model (Outputs Only)');

figure;
semilogx(u, xss,  '-');
legend(model.xNames);
xlabel('Ligand Concentration (L)');
ylabel('Species Concentration');
title('Dose Response Curve for One-Step Model (All species)');