% Copyright (c) 2015, Domenico L. Gatti
% All rights reserved.
% 
% Redistribution and use in source and binary forms, with or without 
% modification, are permitted provided that the following conditions are 
% met:
% 
%     * Redistributions of source code must retain the above copyright 
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright 
%       notice, this list of conditions and the following disclaimer in 
%       the documentation and/or other materials provided with the 
%       distribution
%       
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
% IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
%% Pharmacokinetics-Pharmacodynamics Model (competitive inhibitor)

% This example shows how to analyze a pharmacokinetics-pharmacodynamics
% model in which a drug binds to a carrier protein in the plasma and to a
% target enzyme in a peripheral tissue.

%% Background
% Pharmacokinetics is the study of the mechanisms of absorption and
% distribution of a drug, the chemical changes of the
% substance in the body, and the effects and routes of excretion of the
% metabolites of the drug. Pharmacokinetics is often studied in conjunction
% with pharmacodynamics, the study of a drug's pharmacological effect on
% the body. When both studies are combined to create a model of the drug
% action we refer to this model as a Pharmacokinetics-Pharmacodynamics
% model, or PKPD model.

%% Aims
%
% * Analyze the features of a PKPD model of a drug that inhibits an enzyme 
% in a peripheral tissue.
%
% * Learn how variations in the parameters that describe the pharmacokinetic 
% model affect the drug action.
%
% * Perform a sensitivity analysis to determine which and to what extent
% model parameters affect the drug action.
%
% * Lear how to optimize the drug dosing in order to obtain a stable effect.
%
%%
% We start by loading the project and adding to the path some useful
% programs:

addpath(genpath('../GENERAL_SCRIPTS_FUNCTIONS'));
sbioloadproject('../TOOLBOXES/PKPD/PKPD');

% or alternatively:

% m1 = sbmlimport('PKPD.xml'); % NOTE: do not use xml file because arrays
                               % are reordered. It would require changing
                               % the columns in the rest of the script.

% We also set a file name to save our simulation results:
% savefile = 'PKPD.mat';

%%
% The project contains a simplified model, |m1|, of the plasma, the
% interstitial fluid, the surrenal gland, and all the other tissues. We can
% also load the model interactively by starting SimBiology desktop with
% |simbiology| and opening the project file PKPD.sbproj. The model
% describes a state of hyperaldosteronism due to a surrenal gland tumor, or
% to intrinsic iperactivity of aldosterone synthetase, leading to a higher
% than normal plasma concentration and its correction by means of a 
% competitive inhibitor of this enzyme.

%%
% <<../PKPD.jpg>>

%%
% We can compare this model with a general model of competitive
% inhibition.

%%
% <<../Competitive_inhibition.jpg>>

%%
% As usual, we first get some information about the model.
sbioselect(m1,'Type','compartment')
sbioselect(m1,'Type','species')
sbioselect(m1,'Type','parameter')
sbioselect(m1,'Type','reaction')

% Let's also take a look at the equations that describe the model:
eqs = getequations(m1)

%%
% We will extract key parameters  and species from the model and save
% them as variables in the Workspace.

PL = sbioselect(m1,'Name','Plasma');
SG = sbioselect(m1,'Name','Surrenal_Gland');
IF = sbioselect(m1,'Name','Int_Fluid');
OT = sbioselect(m1,'Name','Other_Tissues');
E = sbioselect(m1,'Name','E');
ES = sbioselect(m1,'Name','ES');    
EI = sbioselect(m1,'Name','EI');
ESI = sbioselect(m1,'Name','ESI');    
Drug_PL = sbioselect(m1,'Name','Drug_Plasma');
Drug_IF = sbioselect(m1,'Name','Drug_Int_fuid');
Drug_SG = sbioselect(m1,'Name','Drug_SG');
Drug_OT = sbioselect(m1,'Name','Drug_Other_Tissues');
Aldosterone_PL = sbioselect(m1,'Name','Aldosterone_Plasma');
Aldosterone_IF = sbioselect(m1,'Name','Aldosterone_IF');
Aldosterone_SG = sbioselect(m1,'Name','Aldosterone_SG');
Aldosterone_OT = sbioselect(m1,'Name','Aldosterone_OT');
Corticosterone_SG = sbioselect(m1,'Name','Corticosterone_SG');
kon_CI = sbioselect(m1,'Name','kon_CI');
koff_CI = sbioselect(m1,'Name','koff_CI');
kon_EI = sbioselect(m1,'Name','kon_EI');
koff_EI = sbioselect(m1,'Name','koff_EI');
kon_E = sbioselect(m1,'Name','kon_E');
koff_E = sbioselect(m1,'Name','koff_E');
kcat_E = sbioselect(m1,'Name','kcat_E');
ksin_cort = sbioselect(m1,'Name','ksin_cort');
kel_cort = sbioselect(m1,'Name','kel_cort');
ksin_E = sbioselect(m1,'Name','ksin_E');
kdeg_E = sbioselect(m1,'Name','kdeg_E');
kdeg_EI = sbioselect(m1,'Name','kdeg_EI');      % degradation EI complex
kdeg_ES = sbioselect(m1,'Name','kdeg_ES');      % degradation ES complex
kti_ald = sbioselect(m1,'Name','kti_ald');
kit_ald = sbioselect(m1,'Name','kit_ald');

% We also store some initial values:
Aldosterone_PL_Init = Aldosterone_PL.InitialAmount;
Aldosterone_SG_Init = Aldosterone_SG.InitialAmount;
Aldosterone_OT_Init = Aldosterone_OT.InitialAmount;
Aldosterone_IF_Init = Aldosterone_IF.InitialAmount;

% We set the koff for the competitive inhibitor to 0.001 (1/hour).
koff_EI.Value = 0.001;

% Kd for the competitive inhibitor.
Drug_Kd = koff_EI.Value/kon_EI.Value;
display(['Drug_Kd = ' num2str(Drug_Kd) ' micrograms/L']);

%% Time evolution of the system
% We will simulate the time evolution of this system. 

% Before we start the simulation we will set all the necessary parameters
% and variable values. We also need to get some information on the
% configuration parameters of the simulation.

cs = getconfigset(m1, 'default');
cs

%% 
% Here we configure the ODE solver:
% cs = getconfigset(m1);
% set(cs, 'SolverType', 'sundials');
% % set(cs, 'SolverType', 'ode15s');
% set(cs, 'StopTime', 6);
% set(cs, 'TimeUnit','hour');
% set(cs.SolverOptions, 'AbsoluteTolerance', 1.e-20);
% set(cs.SolverOptions, 'RelativeTolerance', 1.e-6);
% set(cs.SolverOptions, 'AbsoluteToleranceScaling', true);
% set(cs.SolverOptions, 'SensitivityAnalysis', false);
% set(cs.RunTimeOptions, 'StatesToLog', 'All');
% set(cs.CompileOptions, 'UnitConversion', true);
% set(cs.CompileOptions, 'DimensionalAnalysis', true);
% cs
%%
% We will first set the initial stop time at 6 hours.
Stop = 6;
set(cs, 'StopTime', Stop);

% In the standard simulation we don't want any drug dosing (first []) and 
% no variations in parameters (second []).

PKPD_Kinetics = sbiosimulate(m1, cs, [], []);
PKPD_Kinetics.DataNames
m1.species

% Species are logged in the following order:

%    1         Plasma            Carrier                 microgram/liter
%    2         Plasma            Drug_Carrier            microgram/liter
%    3         Plasma            Drug_Plasma             microgram/liter
%    4         Plasma            Aldosterone_Plasma      microgram/liter
%    5         Int_Fluid         Drug_Int_Fluid          microgram/liter
%    6         Int_Fluid         Aldosterone_IF          microgram/liter
%    7         Surrenal_Gland    E                       microgram/liter
%    8         Surrenal_Gland    Drug_SG                 microgram/liter
%    9         Surrenal_Gland    EI                      microgram/liter
%    10        Surrenal_Gland    Corticosterone_SG       microgram/liter
%    11        Surrenal_Gland    Aldosterone_SG          microgram/liter
%    12        Surrenal_Gland    ES                      microgram/liter
%    13        Surrenal_Gland    ESI                     microgram/liter
%    14        Other_Tissues     Drug_Other_Tissues      microgram/liter
%    15        Other_Tissues     Aldosterone_OT          microgram/liter
   
% We will first be plotting the aldosterone concentration in the plasma and
% the interstitial fluid, and in the tissues where aldosterone is not
% produced.

%%
FIG_1 = figure;
    set(FIG_1,'Units','normalized','Position',[0.6 0.6 0.4 0.4],...
        'Name','PKPD without drug');clf
    
% axes1 = axes('Parent',FIG_1,...
%     'Position',[0.08 0.1 0.9 0.82],'YScale','log',...
%         'XMinorTick','on','XMinorGrid','on');
axes1 = axes('Parent',FIG_1,...
    'Position',[0.08 0.1 0.9 0.82]);

hold(axes1,'all');
box(axes1,'on');
plot(PKPD_Kinetics.Time,PKPD_Kinetics.Data(:,[4 6 11 15]))
legend('Aldosterone\_Plasma','Aldosterone\_IF',...
    'Aldosterone\_SG','Aldosterone\_OT','Location','Best')
% set(gca,'YLim',[-0.05 1.05]);
% set(gca,'XLim',[-0.00001 1.05*Stop]);
ylabel('Species Concentration (micrograms/liter)')
xlabel('Time (hours)')
title('Progress Curve without drug');

%%
% 6 hours is the minimum amount of time required for the system
% to reach equilibrium starting from concentrations of aldosterone that are
% 0 everywhere. Notice how the concentration of aldosterone is rapidly
% equilibrating between plasma and interstitial fluid (IF), but more slowly
% between plasma/(IF) and the tissues. As expected, aldosterone concentration
% is higher in the surrenal gland, where it's synthesized, than everywhere
% else.

%% Sensitivity Analysis - MCA and Control Coefficients 
% So far we have been interested in the dynamic behavior of the model with
% respect to the plasma comcentration of aldosterone. A natural question to
% ask is, which parameters of the model affect aldosterone levels, and what
% are the magnitudes of those effects? To answer this question we can
% compute the sensitivity of plasma aldosterone with respect to various
% parameters in the model. Calculating sensitivities lets you determine
% which species or parameter in a model is most sensitive to a specific
% condition (for example, a drug), defined by a species or parameter. Thus,
% if a model has a species x, and two parameters y and z, the
% time-dependent sensitivities of x with respect to each parameter value
% are the time-dependent derivatives where, the numerator is the
% sensitivity output and the denominators are the sensitivity inputs to
% sensitivity analysis. If we are calculating the sensitivity of plasma
% aldosterone (x)  with respect to two parameters (y) and (z) in the model,
% the time-dependent derivatives are:

%    dx/dy, dx/dz

%%
% Finally, normalizing the sensitivities of a species with respect to
% different parameters allows them to be compared with each other. The
% following lines show you how sensitivities of a species x with respect to
% a parameter k are calculated for each normalization type:
% 
% * 'None' specifies no normalization.
%
%   dx(t)/dk
% 
% * 'Half' specifies normalization relative to the numerator (species
% quantity) only.
%
%   (1/x(t)(dx(t)/dk)  Percentage change of x with respect to a change of k
% 
% * 'Full' specifies that the data should be made dimensionless.
%
%   (k/x(t)(dx(t)/dk)  Percentage change of x with respect to a percentage change of k
%
%
% Related to the concept of sensitivity is the concept of control
% coefficient (cc). A control coefficient is a relative measure of how much a
% perturbation on a 'local' variable affects a 'system' variable (e.g.
% fluxes or concentrations). It is defined as:
% 
% dA/dv(i) * v(i)/A
%

%%
% where A is the variable, i the step (enzyme) and v the rate of the step
% to be perturbed. The most common control coefficients are those for the
% entire pathway flux and the species concentrations, but any variable of
% the system can be analyzed with MCA and have control coefficients defined
% by equations analogous to the above. Although there is no need for the
% system to be in a steady state, traditionally we are interested in
% knowing how the flux J through an entire or part of a metabolic pathway
% (this is the 'system' variable) reacts to a small perturbation of the
% activity of a certain enzyme (this is the 'local' variable) in the
% pathway. 
%
% We will start by calculating the flux (J) control coefficient for the
% model parameters. For example we are interested in calculating the flux J
% of aldosterone synthesis in the surrenal gland; the flux J is the product
% of the time derivative of aldosterone concentration and the surrenal
% gland volume:
%
% J_ald = ReactionFlux19 = (kcat_E*ES)*Surrenal_Gland_Volume
% We also recall that:
% ES is column 12 in array Data
%
% The calculation for the concentration control coefficients is done in the
% same way.

%% 
% We can easily set up separate determinations of flux and concentration
% cc's for all parameters with a single loop:

% Configure task specific stop time.
set(cs, 'StopTime', 30.0); % We simulate up to 30 hours, but we may use less.

params = sbioselect(m1,'Type','parameter');
reactions = sbioselect(m1,'Type','reaction');
species = sbioselect(m1,'Type','species');

% We check parameters and species:
params
PKPD_sens = sbiosimulate(m1);
PKPD_sens.DataNames
m1.species

%%
nparams = size(params,1);
cc_Jk = zeros(1,nparams);
cc_Ck = zeros(1,nparams);

% We can try to make the change small enough to reduce the truncation
% error without increasing the difference error
delta = eps^(1/4);
    
for n = 1:nparams 
    param_no = n; 
    k_ref = params(param_no).Value;
    % Centered difference
    k_range = [(1-delta)*k_ref k_ref (1+delta)*k_ref];
    % 2nd forward difference
    % k_range = [k_ref (1+delta)*k_ref (1+2*delta)*k_ref];
    J_range = zeros(1,3);
    C_range = zeros(1,3);

    for i = 1:3
        params(param_no).Value = k_range(i);
        PKPD_sens = sbiosimulate(m1);
        ReactionFlux19 = (kcat_E.Value*PKPD_sens.Data(end,12))*SG.capacity;
        J_range(i) = ReactionFlux19;                % Flux cc
        C_range(i) = PKPD_sens.Data(end,4);         % Concentration cc
    end

%--------------------------------------------------------------------------    
%     % For dJdv we need an additional calculation because we also need to
%     % record the change in local flux for a small change of k at the same
%     % concentration of the species (involved with k) which is present at 
%     % the original equilibrium. Notice that in this case 'PKPD_sens' is
%     % calculated only once outside the loop.
%     params(param_no).Value = k_ref;
%     PKPD_sens = sbiosimulate(m1);
%     for i = 1:3
%         params(param_no).Value = k_range(i);
%         % The following could be a more complex function involving other
%         % parameters that remain constant.
%         v_range(i) = params(param_no).Value*...
%             PKPD_sens.Data(end,'variable affected by param_no');
%     end
%     % Centered difference
%     dJdv = ((J_range(3)-J_range(1))/(v_range(3)-v_range(1)))*(v_range(2)/J_range(2)); 
%     cc_Jv(n) = dJdv;
%     % If the flux is a function of only k, then there is no difference in
%     % the final result with respect to dJdk, but if the flux is a more
%     % complex function involving other parameters, then dJdv can be vary
%     % different from dJdk.
%--------------------------------------------------------------------------    
    
% 2nd forward difference
%     dJdk = ((-J_range(3)+4*J_range(2)-3*J_range(1))/(k_range(3)-k_range(1))) * k_range(1)/J_range(1); 
%     dCdk = ((-C_range(3)+4*C_range(2)-3*C_range(1))/(k_range(3)-k_range(1))) * k_range(1)/C_range(1);
% Centered difference
    dJdk = ((J_range(3)-J_range(1))/(k_range(3)-k_range(1)))*(k_range(2)/J_range(2)); 
    dCdk = ((C_range(3)-C_range(1))/(k_range(3)-k_range(1)))*(k_range(2)/C_range(2)); 
    % end
    cc_Jk(n) = dJdk;  
    cc_Ck(n) = dCdk;  
    params(param_no).Value = k_ref;
end

%%

CONTROL_COEFF = figure;
    set(CONTROL_COEFF,'Units','normalized','Position',[0.4 0.2 0.8 0.8],...
        'Name','Control Coefficients');clf

% Flux control coefficients
nz_ind = (find(cc_Jk));
nnames = length(nz_ind);
names = cell(nnames,1);
for i = 1:nnames
    names{i} = params(nz_ind(i)).Name;
end
    
    haxes1 = axes('Parent',CONTROL_COEFF,...
    'Position',[0.058 0.53 0.43 0.42]);

        hbar = barh((cc_Jk(nz_ind)));
%        haxes1 = get(hbar(1), 'Parent');
%         set(haxes,'XScale','log',...
%         'XMinorTick','on','XMinorGrid','on');
        set(haxes1, 'ytick', 1:nnames);
        set(haxes1, 'yticklabel', names);
        set(haxes1, 'XGrid', 'on');
        xlabel('cc(J\_ald,parameter)');
        ylabel('Parameters');
        title('Flux Control Coefficients');

% Zoom in for weak signals        
clear nz_ind nnames names        
nz_ind = (find(abs(cc_Jk) > 0 & abs(cc_Jk) <= 0.1));
nnames = length(nz_ind);
names = cell(nnames,1);
for i = 1:nnames
    names{i} = params(nz_ind(i)).Name;
end
               
    haxes2 = axes('Parent',CONTROL_COEFF,...
    'Position',[0.058 0.06 0.43 0.42]);
        
        hbar = barh((cc_Jk(nz_ind)));
%        haxes2 = get(hbar(1), 'Parent');
%         set(haxes,'XScale','log',...
%         'XMinorTick','on','XMinorGrid','on');
        set(haxes2, 'ytick', 1:nnames);
        set(haxes2, 'yticklabel', names);
        set(haxes2, 'XGrid', 'on');
        xlabel('cc(J\_ald,parameter)');
        ylabel('Parameters');
        % set(gca,'XLim',[-0.002 0.002]);

% Concentration control coefficients                
clear nz_ind nnames names        
nz_ind = (find(cc_Ck));
nnames = length(nz_ind);
names = cell(nnames,1);
for i = 1:nnames
    names{i} = params(nz_ind(i)).Name;
end
        
    haxes3 = axes('Parent',CONTROL_COEFF,...
    'Position',[0.55 0.53 0.43 0.42]);

        hbar = barh((cc_Ck(nz_ind)));
%        haxes1 = get(hbar(1), 'Parent');
%         set(haxes,'XScale','log',...
%         'XMinorTick','on','XMinorGrid','on');
        set(haxes3, 'ytick', 1:nnames);
        set(haxes3, 'yticklabel', names);
        set(haxes3, 'XGrid', 'on');
        xlabel('cc(PL\_ald,parameter)');
        ylabel('Parameters');
        % set(gca,'XLim',[-1.1 1.1]);
        title('Concentration Control Coefficients');

% Zoom in for weak signals        
clear nz_ind nnames names        
nz_ind = (find(abs(cc_Ck) > 0 & abs(cc_Ck) <= 0.1));
nnames = length(nz_ind);
names = cell(nnames,1);
for i = 1:nnames
    names{i} = params(nz_ind(i)).Name;
end
                       
    haxes4 = axes('Parent',CONTROL_COEFF,...
    'Position',[0.55 0.06 0.43 0.42]);
        
        hbar = barh((cc_Ck(nz_ind)));
%        haxes2 = get(hbar(1), 'Parent');
%         set(haxes,'XScale','log',...
%         'XMinorTick','on','XMinorGrid','on');
        set(haxes4, 'ytick', 1:nnames);
        set(haxes4, 'yticklabel', names);
        set(haxes4, 'XGrid', 'on');
        xlabel('cc(PL\_ald,parameter)');
        ylabel('Parameters');
        % set(gca,'XLim',[-0.002 0.002]);

%%
% Negative control coefficients indicate that an increase in that parameter
% produces a decrease in flux.

%% 
% % We can also calculate the ccJ with respect to combinations of variables
% % perturbed. For example we could calculate the flux control coefficient on
% % the sythesis of aldosterone (cc_Jald) of kon and koff acting together and
% % leaving the Kd for corticosterone unchanged:
% 
% param1_no = 18; % This is kon_E
% param2_no = 19; % This is koff_E
% delta = 0.0001;
% k_ref1 = params(param1_no).Value;
% k_ref2 = params(param2_no).Value;
% k_range1 = [(1-delta)*k_ref1 k_ref1 (1+delta)*k_ref1];
% k_range2 = [(1-delta)*k_ref2 k_ref2 (1+delta)*k_ref2];
% 
% J_range = zeros(1,3);
% v_range = zeros(1,3);
% 
% % First the J_range: PKPD_sens is inside the loop.
% for i = 1:3
%     params(param1_no).Value = k_range1(i);
%     params(param2_no).Value = k_range2(i);
%     PKPD_sens = sbiosimulate(m1);
%     ReactionFlux19 = (kcat_E.Value*PKPD_sens.Data(end,12))*SG.capacity;
%     J_range(i) = ReactionFlux19;
% end
% 
% % Next the v_range: PKPD_sens is outside the loop.
%     params(param1_no).Value = k_ref1;
%     params(param2_no).Value = k_ref2;
%     PKPD_sens = sbiosimulate(m1);
%     for i = 1:3
%         params(param1_no).Value = k_range1(i);
%         params(param2_no).Value = k_range2(i);
%         % The following is the flux of the reaction of corticosterone
%         % binding to the E.
%         v_range(i) = params(param1_no).Value*PKPD_sens.Data(end,10)*...
%             PKPD_sens.Data(end,7)*SG.Capacity - ...
%             params(param2_no).Value*PKPD_sens.Data(end,12)*SG.Capacity;
%     end
% 
%     dJdk1 = ((J_range(3)-J_range(1))/(k_range1(3)-k_range1(1)))*(k_range1(2)/J_range(2)); 
%     dJdk2 = ((J_range(3)-J_range(1))/(k_range2(3)-k_range2(1)))*(k_range2(2)/J_range(2)); 
%     dJdv = ((J_range(3)-J_range(1))/(v_range(3)-v_range(1)))*(v_range(2)/J_range(2)); 
% 
% 
% % Before we forget we return kon_AS to its original value:
% params(param1_no).Value = k_ref1;
% params(param2_no).Value = k_ref2;
% 
% % 
% ccJk1 = dJdk1
% ccJk2 = dJdk2
% ccJv12 = dJdv
% 
% 
% % We can see how the ccJ of either kon_E or koff_E when both are changed
% % simultaneously is much smaller than the individual ccJ's. As a general
% % rule the combined effect of two ccJ on a step is approximately equal to
% % the sum of the individual ccJ's. In this case, since they have opposite
% % sign, they tend to almost complitely obliterate each other.
% %
%% Dose Response Curves
% It is clear that the rate of synthesis of corticosterone (Cort syn) and
% the rate of renal clearance of aldosterone (Ald el) are the main factors
% affecting the plasma concemtration of aldosterone. Unfortunately we
% cannot increase the renal clearance and we do not have a drug that
% inhibits corticosterone synthesis. The third most important factor is the
% kcat (kcat_E) of aldosterone synthase, but we only have a competitive
% inhibitor that does not affect the kcat, and in practive it works by
% subtracting a certain amount of enzyme from the pool of active synthase.
% We will work with what we have.

% The first step is to determine what level of inhibition we need to
% achieve in order to bring the plasma aldosterone level around 3
% micrograms/liter. We can set up a Dose-Response curve to determine this
% level.

E_ref = E.InitialAmount;

High_Conc = 1000;
ncons = 35;
conc_vec = logspace(log10(1),log10(High_Conc),ncons);
% openvar conc_vec
range = conc_vec <= 1 + High_Conc;

% We simulate for 30 hours
Stop = 30;
set(cs, 'StopTime', Stop);
Response = zeros(1,ncons);
Dose = conc_vec;

for i = 1:ncons
    E.InitialAmount = conc_vec(i);
    PKPD_Kinetics = sbiosimulate(m1, cs, [], []);
    Response(i) = PKPD_Kinetics.Data(end,4);
    % Dose(i) = PKPD_Kinetics.Data(end,7);
end

% Before we forget, we return the Aldosterone Synthase concentration to the
% original value:
E.InitialAmount = E_ref;

%%
E_DOSE_RESPONSE = figure;
    set(E_DOSE_RESPONSE,'Units','normalized','Position',[0.6 0.6 0.4 0.4],...
        'Name','Dose Response Curve');clf
axes_edr = axes('Parent',E_DOSE_RESPONSE,...
    'Position',[0.08 0.1 0.86 0.82],'XScale','log',...
        'XMinorTick','on','XMinorGrid','on','XDir','reverse');
hold(axes_edr,'all');
box(axes_edr,'on');

% Dose-Response plots (percentage of maximal response against log of the
% drug concentration):

semilogx(Dose,Response,'-or',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','y',...
             'MarkerSize',5,'Parent',axes_edr);

target_E = Dose(find(Response == min(abs(Response-0.3)+0.3))); 

vline([target_E target_E],{'g-', 'LineWidth', 3});
% hline([ly ly],'g--');         
% set(gca,'YLim',[-0.05 1.1]);
% set(gca,'XLim',[-5 5+conc_vec(end)]);
% set(gca,'XLim',[-5 50]);
         
legend('Plasma Ald','Location','Best');
xlabel('Aldosterone Synthase (micrograms/L)');
ylabel('Plasma aldosterone at 30 hour (micrograms/L)');
title('Dose-Response Curve for AS concentration');

string1 = 'Target [E] = ';
string2 = num2str(target_E,'%6.3f\n');
string3 = ' microgram/L';

% Create textbox
annotation(E_DOSE_RESPONSE,'textbox',...
    [0.55 0.72 0.35 0.100],...
    'String',{[string1 string2 string3]},...
    'FontWeight','bold',...
    'FontSize',14,...
    'FitBoxToText','off',...
    'BackgroundColor',[1 0.968627452850342 0.921568632125854],...
    'Color',[1 0 0]);

%%
% We are almost there: this means that we need to bring the activity of
% aldosterone synthase to ~4.7% of the initial amount in order to
% bring the concentration of plasma aldosterone into ~normal range. Now we
% are ready to start experimenting with our drug.

% We get some information about the dosing that is already stored in 
% the model:
doses = get(m1.doses);
d1 = getdose(m1, 'Single Dose');
d2 = getdose(m1, 'Repeated Dose');

% Next we develop a Dose-Response curve of the level of plasma aldosterone
% at 30 hours after a single administration at 6 hour. As usual we will
% loop through a vector of concentrations spaced logarithmically.

High_Conc = 2500;
ncons = 35;
conc_vec = logspace(log10(1),log10(High_Conc),ncons);
% openvar conc_vec
range = conc_vec <= 1 + High_Conc;

% We simulate for 30 hours
Stop = 30;
set(cs, 'StopTime', Stop);
Response = zeros(1,ncons);
Dose = conc_vec;
Drug_conc = zeros(1,ncons);

for i = 1:ncons
    set(d1, 'Amount', conc_vec(i));
    set(d1, 'Rate', conc_vec(i));
    PKPD_Kinetics = sbiosimulate(m1, cs, [], d1);
    Response(i) = PKPD_Kinetics.Data(end,4);
    Drug_conc(i) = PKPD_Kinetics.Data(end,8);
end

%%
DOSE_RESPONSE = figure;
    set(DOSE_RESPONSE,'Units','normalized','Position',[0.6 0.2 0.4 0.8],...
        'Name','Dose Response Curve');clf
axes_dr1 = axes('Parent',DOSE_RESPONSE,...
    'Position',[0.08 0.53 0.86 0.42],'XScale','log',...
        'XMinorTick','on','XMinorGrid','on');
hold(axes_dr1,'all');
box(axes_dr1,'on');

% Dose-Response plots (percentage of maximal response against log of the
% drug concentration):

semilogx(Dose,Response,'or',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','y',...
             'MarkerSize',5,'Parent',axes_dr1);

% hline([ly ly],'g--');         
% set(gca,'YLim',[-0.05 1.1]);
% set(gca,'XLim',[-5 5+conc_vec(end)]);
% set(gca,'XLim',[-5 50]);

%--------------------------------------------------------------------------
% A Dose-Response Logistic fit is based on the following equation:
% Response = (Max-Min)/(1+(Kd/x)^n  in which n = nH = Hill coefficient. 
% Notice that theoretically nH is not supposed to be any higher than the
% number of binding sites which are occupied simultaneously. In the
% equation below a = ID50 and b = nH.
%--------------------------------------------------------------------------

Y = Response;
X = Dose;

% We use a 4 parameter logistic fit:
         
f = fittype('(c - d)/(1 + (a/x)^b)');
[Hill_1,GOF_1] = fit(X',Y',f,'StartPoint',[ 80 -1 0.8 0.05 ]);
Hill_1

semilogx(Dose,Hill_1(Dose)','-r',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','y',...
             'MarkerSize',5,'Parent',axes_dr1);
         
legend('Plasma Ald','Logistic Fit','Location','Best');
xlabel('Drug Dose (mg) at 6 hour');
ylabel('Plasma aldosterone at 30 hour (micrograms/l)');
title('Dose-Response Curve for AS inhibition');

fit_params = coeffvalues(Hill_1);
ID50 = fit_params(1);
string1 = 'ID50 = ';
string2 = num2str(ID50,'%6.3f\n');
string3 = ' mg';

% Create textbox
annotation(DOSE_RESPONSE,'textbox',...
    [0.70 0.82 0.18 0.06],...
    'String',{[string1 string2 string3]},...
    'FontWeight','bold',...
    'FontSize',14,...
    'FitBoxToText','off',...
    'BackgroundColor',[1 0.968627452850342 0.921568632125854],...
    'Color',[1 0 0]);



axes_dr2 = axes('Parent',DOSE_RESPONSE,...
    'Position',[0.08 0.06 0.86 0.42],'XScale','log',...
        'XMinorTick','on','XMinorGrid','on');
hold(axes_dr2,'all');
box(axes_dr2,'on');

% Dose-Response plots (percentage of maximal response against log of the
% ligand concentration):

semilogx(Drug_conc,Response,'ob',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','y',...
             'MarkerSize',5,'Parent',axes_dr2);

% set(gca,'YLim',[-0.05 1.1]);
% set(gca,'XLim',[-5 5+conc_vec(end)]);
% set(gca,'XLim',[-5 50]);

%--------------------------------------------------------------------------
% Logistic fit In the
% equation below a = Ki(app) and b = nH.

Y = Response;
X = Drug_conc;

% We use a 4 parameter logistic fit:
         
f = fittype('(c - d)/(1 + (a/x)^b)');
[Hill_2,GOF_2] = fit(X',Y',f,'StartPoint',[ 1800 -1 0.82 0.03 ],...
    'Lower',[1000 -1.4 0.72 0.01],'Upper',[3000 -0.6 0.92 0.05]);
Hill_2

% plot(Hill_2,'--b','Parent',axes_dr2);

semilogx(Drug_conc,Hill_2(Drug_conc)','-b',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','y',...
             'MarkerSize',5,'Parent',axes_dr2);

legend('Plasma Ald','Logistic Fit','Location','Best');
xlabel('Surrenal Drug Concentration (micrograms/l) at 30 hour');
ylabel('Plasma aldosterone at 30 hour (micrograms/l)');
% title('Dose-Response Curve for AS inhibition');

fit_params = coeffvalues(Hill_2);
Ki_app = fit_params(1);
string1 = 'Ki\_app = ';
string2 = num2str(Ki_app,'%6.3f\n');
string3 = ' micrograms/l';

% Create textbox
annotation(DOSE_RESPONSE,'textbox',...
    [0.58 0.34 0.34 0.06],...
    'String',{[string1 string2 string3]},...
    'FontWeight','bold',...
    'FontSize',14,...
    'FitBoxToText','off',...
    'BackgroundColor',[1 0.968627452850342 0.921568632125854],...
    'Color',[0 0 1]);

%%
% The dose response curves suggest a ID50 (dose inhibiting 50% of the
% enzyme) around 80 mg, and an apparent Ki for the drug of 1.4 mg/l. This is
% much higher than the expected Ki = Kd = 0.01 micrograms/l. We will return to
% this unusual result in the example 'PKPDnc'.

%% Developing a Drug Dosing Schedule
% Based on the ID50(at 30 hours) of 80 mg and in consideration of the
% progressive accumulation of the drug with repeated doses, we can try a
% schedule of 30 mg/day. Let's see what the stored repeated dosage looks
% like:

get(d2)

%%
% Let's set the repeat count to 20 and the amount to 30 mg:

d2.RepeatCount = 20;
d2.Amount = 30;
d2.Rate = 30;

% It looks reasonable: now we can repeat the simulation with this dose
% schedule and simulate for 960 hours (40 days). After 20 days we will see
% the aldosterone level going up again, but we are primarily interested in
% seeing if during the administration of the drug we obtain a stable normal
% level of plasma aldosterone

Stop = 960;
set(cs, 'StopTime', Stop);

PKPD_Kinetics = sbiosimulate(m1, cs, [], d2);

%%
FIG_6 = figure;
    set(FIG_6,'Units','normalized','Position',[0.6 0.6 0.4 0.4],...
        'Name','PKPD with drug');clf
    
% axes1 = axes('Parent',FIG_1,...
%     'Position',[0.08 0.1 0.9 0.82],'YScale','log',...
%         'XMinorTick','on','XMinorGrid','on');
axes6 = axes('Parent',FIG_6,...
    'Position',[0.08 0.1 0.9 0.82]);

hold(axes6,'all');
box(axes6,'on');
plot(PKPD_Kinetics.Time,PKPD_Kinetics.Data(:,[4 6 15]))
legend('Aldosterone\_Plasma','Aldosterone\_IF','Aldosterone\_OT',...
    'Location','Best')
% set(gca,'YLim',[-0.05 1.05]);
% set(gca,'XLim',[-0.00001 1.05*Stop]);
ylabel('Species Concentration (microgram/liter)')
xlabel('Time (hours)')
title('Progress Curve with 30 mg/day drug');

%%
% That looks very good. We can probably achieve an optimal result with a
% repeated larger dose taken 1/day. 
d2.RepeatCount = 20;
d2.Amount = 60;
d2.Interval = 24;
d2.Rate = 60;
d2.StartTime = 6;

% It looks reasonable: now we can repeat the simulation with this dose
% schedule and simulate for 960 hours (40 days). After 20 days we will see
% the aldosterone level going up again, but we are primarily interested in
% seeing if during the administration of the drug we obtain a stable normal
% level of plasma aldosterone

Stop = 960;
set(cs, 'StopTime', Stop);

PKPD_Kinetics = sbiosimulate(m1, cs, [], d2);

%%
FIG_7 = figure;
    set(FIG_7,'Units','normalized','Position',[0.6 0.6 0.4 0.4],...
        'Name','PKPD with drug');clf
    
% axes1 = axes('Parent',FIG_1,...
%     'Position',[0.08 0.1 0.9 0.82],'YScale','log',...
%         'XMinorTick','on','XMinorGrid','on');
axes7 = axes('Parent',FIG_7,...
    'Position',[0.08 0.1 0.9 0.82]);

hold(axes7,'all');
box(axes7,'on');
plot(PKPD_Kinetics.Time,PKPD_Kinetics.Data(:,[4 6 15]))
legend('Aldosterone\_Plasma','Aldosterone\_IF','Aldosterone\_OT',...
    'Location','Best')
% set(gca,'YLim',[-0.05 1.05]);
% set(gca,'XLim',[-0.00001 1.05*Stop]);
ylabel('Species Concentration (micrograms/liter)')
xlabel('Time (hours)')
title('Progress Curve with 60 mg/day for 20 days ');

%% Conclusions
% This looks fantastic! We have achieved a stable normalization of the
% plasma aldosterone levels with only 1 dayly dose and acceptably small 
% daily fluctuations. We can expect that even smaller fluctuations could be
% achieved by splitting the daily dose in two doses to be taken every 12
% hours.

% Let's decrease the amount to 25 mg and increase the repeat count to 40 
% with an interval of 12 hours:

d2.Amount = 25;
d2.Rate = 25;
d2.Interval = 12;
d2.RepeatCount = 40;
d2.StartTime = 6;

%%

Stop = 960;
set(cs, 'StopTime', Stop);

PKPD_Kinetics = sbiosimulate(m1, cs, [], d2);

%%
FIG_8 = figure;
    set(FIG_8,'Units','normalized','Position',[0.6 0.6 0.4 0.4],...
        'Name','PKPD with drug');clf
    
% axes1 = axes('Parent',FIG_1,...
%     'Position',[0.08 0.1 0.9 0.82],'YScale','log',...
%         'XMinorTick','on','XMinorGrid','on');
axes8 = axes('Parent',FIG_8,...
    'Position',[0.08 0.1 0.9 0.82]);

hold(axes8,'all');
box(axes8,'on');
plot(PKPD_Kinetics.Time,PKPD_Kinetics.Data(:,[4 6 15]))
hline([0.34 0.34],{'g-','LineWidth',2});
legend('Aldosterone\_Plasma','Aldosterone\_IF','Aldosterone\_OT',...
    'Location','Best')
% set(gca,'YLim',[-0.05 1.05]);
% set(gca,'XLim',[-0.00001 1.05*Stop]);
ylabel('Species Concentration (micrograms/liter)')
xlabel('Time (hours)')
title('Progress Curve with 25 mg twice/day for 20 days ');

%% 
% Finally, we are interested in knowing whether the cc's for the
% aldosterone synthesis flux have changed in this new steady state produced
% by the drug. Again we can easily determine the cc's for all parameters
% with a single loop. A convenient time to stop the simulation is 398 hours
% (middle of a small fluctuation):

Stop = 397;
set(cs, 'StopTime', Stop);

nparams = size(params,1);
cc_Jk = zeros(1,nparams);
cc_Ck = zeros(1,nparams);
delta = eps^(1/4);

for n = 1:nparams 
    param_no = n; 
    k_ref = params(param_no).Value;
    k_range = [(1-delta)*k_ref k_ref (1+delta)*k_ref];
    J_range = zeros(1,3);
    C_range = zeros(1,3);

    for i = 1:3
        params(param_no).Value = k_range(i);
        PKPD_sens = sbiosimulate(m1, cs, [], d2);
        ReactionFlux19 = (kcat_E.Value*PKPD_sens.Data(end,12))*SG.capacity;
        J_range(i) = ReactionFlux19;                % Flux cc
        C_range(i) = PKPD_sens.Data(end,4);         % Concentration cc
    end
    
    dJdk = ((J_range(3)-J_range(1))/(k_range(3)-k_range(1)))*(k_range(2)/J_range(2)); 
    dCdk = ((C_range(3)-C_range(1))/(k_range(3)-k_range(1)))*(k_range(2)/C_range(2)); 
    cc_Jk(n) = dJdk;  
    cc_Ck(n) = dCdk;      
    params(param_no).Value = k_ref;
end

%%
% and derive the value of each cc

CONTROL_COEFF_2 = figure;
    set(CONTROL_COEFF_2,'Units','normalized','Position',[0.4 0.2 0.8 0.8],...
        'Name','Control Coefficients');clf

% Flux control coefficients
nz_ind = (find(cc_Jk));
nnames = length(nz_ind);
names = cell(nnames,1);
for i = 1:nnames
    names{i} = params(nz_ind(i)).Name;
end
    
    haxes1 = axes('Parent',CONTROL_COEFF_2,...
    'Position',[0.058 0.53 0.43 0.42]);

        hbar = barh((cc_Jk(nz_ind)));
%        haxes1 = get(hbar(1), 'Parent');
%         set(haxes,'XScale','log',...
%         'XMinorTick','on','XMinorGrid','on');
        set(haxes1, 'ytick', 1:nnames);
        set(haxes1, 'yticklabel', names);
        set(haxes1, 'XGrid', 'on');
        xlabel('cc(J\_ald,parameter)');
        ylabel('Parameters');
        title('Flux Control Coefficients');

% Zoom in for weak signals        
clear nz_ind nnames names        
nz_ind = (find(abs(cc_Jk) > 0 & abs(cc_Jk) <= 0.1));
nnames = length(nz_ind);
names = cell(nnames,1);
for i = 1:nnames
    names{i} = params(nz_ind(i)).Name;
end
        
    haxes2 = axes('Parent',CONTROL_COEFF_2,...
    'Position',[0.058 0.06 0.43 0.42]);
        
        hbar = barh((cc_Jk(nz_ind)));
%        haxes2 = get(hbar(1), 'Parent');
%         set(haxes,'XScale','log',...
%         'XMinorTick','on','XMinorGrid','on');
        set(haxes2, 'ytick', 1:nnames);
        set(haxes2, 'yticklabel', names);
        set(haxes2, 'XGrid', 'on');
        xlabel('cc(J\_ald,parameter)');
        ylabel('Parameters');
        set(gca,'XLim',[-0.03 0.04]);

% Concentration control coefficients                
clear nz_ind nnames names        
nz_ind = (find(cc_Ck));
nnames = length(nz_ind);
names = cell(nnames,1);
for i = 1:nnames
    names{i} = params(nz_ind(i)).Name;
end
        
    haxes3 = axes('Parent',CONTROL_COEFF_2,...
    'Position',[0.55 0.53 0.43 0.42]);

        hbar = barh((cc_Ck(nz_ind)));
%        haxes1 = get(hbar(1), 'Parent');
%         set(haxes,'XScale','log',...
%         'XMinorTick','on','XMinorGrid','on');
        set(haxes3, 'ytick', 1:nnames);
        set(haxes3, 'yticklabel', names);
        set(haxes3, 'XGrid', 'on');
        xlabel('cc(PL\_ald,parameter)');
        ylabel('Parameters');
        set(gca,'XLim',[-1.1 1.1]);
        title('Concentration Control Coefficients');

% Zoom in for weak signals        
clear nz_ind nnames names        
nz_ind = (find(abs(cc_Ck) > 0 & abs(cc_Ck) <= 0.1));
nnames = length(nz_ind);
names = cell(nnames,1);
for i = 1:nnames
    names{i} = params(nz_ind(i)).Name;
end
        
    haxes4 = axes('Parent',CONTROL_COEFF_2,...
    'Position',[0.55 0.06 0.43 0.42]);
        
        hbar = barh((cc_Ck(nz_ind)));
%        haxes2 = get(hbar(1), 'Parent');
%         set(haxes,'XScale','log',...
%         'XMinorTick','on','XMinorGrid','on');
        set(haxes4, 'ytick', 1:nnames);
        set(haxes4, 'yticklabel', names);
        set(haxes4, 'XGrid', 'on');
        xlabel('cc(PL\_ald,parameter)');
        ylabel('Parameters');
        set(gca,'XLim',[-0.03 0.04]);

%%
% The presence of the inhibitor dramatically affects the ccJ for the
% synthesis of aldosterone. Thus, other steps in the kinetic models become
% susceptible of further modulation.

%%
close all
% save(savefile);
