% 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.
%
%% General dependencies
% We always start from the PRACTICE directory and we add to the path
% subdirectories containing various tools described in the book chapters.
addpath(genpath('../GENERAL_SCRIPTS_FUNCTIONS'));
% addpath(genpath('../DATABASE'));
% addpath(genpath('../TOOLBOXES'));
addpath(genpath('../TOOLBOXES/PKPD'));

%% Pharmacokinetics-Pharmacodynamics Model (uncompetitive drug)

% This example shows how to analyze a pharmacokinetics-pharmacodynamics
% model in which a drug, acting as a uncompetitive inhibitor, binds to a
% carrier protein in the plasma and to a target protein 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 an uncompetitive drug that
% inhibits an enzyme in a peripheral tissue.
%
% * Learn how variations in the parameters that describe the pharmacokinetic 
% model affect the drug action.
%
% * Lear how to optimize the drug dosing in order to obtain a stable effect.
%
%%
sbioloadproject('PKPDunc');

%%
% 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 an
% uncompetitive inhibitor of this enzyme.

%%
% <<../PKPDunc.jpg>>

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

%%
% <<../Uncompetitive_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');
Drug_PL = sbioselect(m1,'Name','Drug_Plasma');
Drug_SG = sbioselect(m1,'Name','Drug_SG');
Drug_IF = sbioselect(m1,'Name','Drug_Int_Fluid');
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');
E = sbioselect(m1,'Name','E');
ES = sbioselect(m1,'Name','ES');    
ESI = sbioselect(m1,'Name','ESI');    
kon = sbioselect(m1,'Name','kon');      % inhibitor binding to carrier
koff = sbioselect(m1,'Name','koff');
kon3 = sbioselect(m1,'Name','kon3');    % inhibitor binding to ES
koff3 = sbioselect(m1,'Name','koff3');
kon_AS = sbioselect(m1,'Name','kon_AS');    % substrate binding to E
koff_AS = sbioselect(m1,'Name','koff_AS');
kcat_AS = sbioselect(m1,'Name','kcat_AS');
ksin_cort = sbioselect(m1,'Name','ksin_cort');
kel_cort = sbioselect(m1,'Name','kel_cort');
ksin2 = sbioselect(m1,'Name','ksin2');
kdeg2 = sbioselect(m1,'Name','kdeg2');
kdeg_ES = sbioselect(m1,'Name','kdeg_ES');
kdeg_ESI = sbioselect(m1,'Name','kdeg_ESI');
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;

Drug_Kd = koff3.Value/kon3.Value
display(['Drug_Kd = ' num2str(Drug_Kd) ' micrograms/L']);

% The Kd for this uncompetitive inhibitor is 100 times bigger than that
% of the competitive inhibitor studied in the previous example.

%% 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

%%
% 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    Aldosterone_Synthase    microgram/liter
%    8         Surrenal_Gland    Drug_SG                 microgram/liter
%    9         Surrenal_Gland    Corticosterone_SG       microgram/liter
%    10        Surrenal_Gland    Aldosterone_SG          microgram/liter
%    11        Surrenal_Gland    ES                      microgram/liter
%    12        Surrenal_Gland    ESI                     microgram/liter
%    13        Other_Tissues     Drug_Other_Tissues      microgram/liter
%    14        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 10 14]))
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.

%% MCA and Control Coefficients 
% We can check the control coefficients:

params = sbioselect(m1,'Type','parameter');

% We check parameters, species, and the names of simulation states logged:
params
PKPD_sens = sbiosimulate(m1);
PKPD_sens.DataNames
m1.species

%% 
% We get separate determinations of cc's for all parameters with a single
% loop:

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);

    for i = 1:3
        params(param_no).Value = k_range(i);
        PKPD_sens = sbiosimulate(m1);
        ReactionFlux17 = (kcat_AS.Value*PKPD_sens.Data(end,11))*SG.capacity;
        J_range(i) = ReactionFlux17;
        C_range(i) = PKPD_sens.Data(end,4);                
    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

%%
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

CONTROL_COEFF_1 = figure;
    set(CONTROL_COEFF_1,'Units','normalized','Position',[0.6 0.6 0.8 0.4],...
        'Name','Integrated Sensitivity Analysis');clf
    subplot(1,2,1)
        hbar = barh((cc_Jk(nz_ind)));
        haxes = get(hbar(1), 'Parent');
%         set(haxes,'XScale','log',...
%         'XMinorTick','on','XMinorGrid','on');
        set(haxes, 'ytick', 1:nnames);
        set(haxes, 'yticklabel', names);
        xlabel('cc(J\_ald,parameter)');
        ylabel('Parameters');
        set(gca,'XLim',[-1.2 1.2]);
        title('Flux Control Coefficients');

    subplot(1,2,2)
        hbar = barh((cc_Ck(nz_ind)));
        haxes = get(hbar(1), 'Parent');
%         set(haxes,'XScale','log',...
%         'XMinorTick','on','XMinorGrid','on');
        set(haxes, 'ytick', 1:nnames);
        set(haxes, 'yticklabel', names);
        xlabel('cc(C\_ald,parameter)');
        ylabel('Parameters');
        set(gca,'XLim',[-1.2 1.2]);
        title('Concentration Control Coefficients');
        
%%    
% Negative control coefficients indicate that an increase in that parameter
% produces a decrease in flux.

%% Changes in the amount of surrenal aldosterone synthase 
% 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 (AS_kcat) of aldosterone synthase.

% 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. Let's see what variations of the model with less enzyme
% are already available.

m1.variants

%    Index:  Name:             Active:
%    1       V1                false
%    2       V2                false
%    3       V3                false
% 
% We will activate the variant A and run a simulation for 6 hours
% (default configuration setting).

% Initialize configset for analysis run.
% cs = getconfigset(m1, 'default');
Stop = 6;
set(cs, 'StopTime', Stop);

% Assign variable name to object.
variant1 = getvariant(m1, 'V1');
set(variant1, 'Content', {{'species', ...
    'Surrenal_Gland.E', 'InitialAmount', 100}});

% Make sure the variant is active
m1.variant(1).Active = 1;
m1.variant(2).Active = 0;
m1.variant(3).Active = 0;

% Initialize variants for analysis run.
variants = variant1;

% RUNSIMULATION simulate SimBiology model, m1 using the
% variant object(s), variants. Variant object settings supercede
% model property-values.

% Run simulation.
PKPD_Kinetics = sbiosimulate(m1, cs, variants, []);

%%
FIG_3 = figure;
    set(FIG_3,'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');
axes3 = axes('Parent',FIG_3,...
    'Position',[0.08 0.1 0.9 0.82]);

hold(axes3,'all');
box(axes3,'on');
plot(PKPD_Kinetics.Time,PKPD_Kinetics.Data(:,[4 6 10 14]))
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 with 1/10 AS');

%%
% Nothing seems to change. We can try 6/100 AS. Since there is
% less enzyme we expect the equilibrium to be reached more slowly; thus we
% will extend the simulation to 10 hours.

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

% Configure properties.
variant2 = getvariant(m1, 'V2');
set(variant2, 'Content', {{'species', ...
    'Surrenal_Gland.E', 'InitialAmount', 60}});

% Make sure the variant is active
m1.variant(1).Active = 0;
m1.variant(2).Active = 1;
m1.variant(3).Active = 0;

% Initialize variants for analysis run.
variants = variant2;

% RUNSIMULATION simulate SimBiology model, m1 using the
% variant object(s), variants. Variant object settings supercede
% model property-values.

% Run simulation.
PKPD_Kinetics = sbiosimulate(m1, cs, variants, []);

%%
FIG_4 = figure;
    set(FIG_4,'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');
axes4 = axes('Parent',FIG_4,...
    'Position',[0.08 0.1 0.9 0.82]);

hold(axes4,'all');
box(axes4,'on');
plot(PKPD_Kinetics.Time,PKPD_Kinetics.Data(:,[4 6 10 14]))
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 with 6/100 AS');

%%
% The equilibrium was reached more slowly, but the final concentration of plasma 
% aldosterone was still too high. We can try 3/100 AS. 
Stop = 10;
set(cs, 'StopTime', Stop);

% Configure properties.
variant3 = getvariant(m1, 'V3');
set(variant3, 'Content', {{'species', ...
    'Surrenal_Gland.E', 'InitialAmount', 30}});

% Make sure the variant is active
m1.variant(1).Active = 0;
m1.variant(2).Active = 0;
m1.variant(3).Active = 1;

% Initialize variants for analysis run.
variants = variant3;

% RUNSIMULATION simulate SimBiology model, m1 using the
% variant object(s), variants. Variant object settings supercede
% model property-values.

% Run simulation.
PKPD_Kinetics = sbiosimulate(m1, cs, variants, []);

%%
FIG_5 = figure;
    set(FIG_5,'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');
axes5 = axes('Parent',FIG_5,...
    'Position',[0.08 0.1 0.9 0.82]);

hold(axes5,'all');
box(axes5,'on');
plot(PKPD_Kinetics.Time,PKPD_Kinetics.Data(:,[4 6 10 14]))
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 with 3/100 AS');

%% Dose Response Curve
% We are almost there: this means that we need to bring the activity of
% aldosterone synthase to ~ 4/100 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.

% First we set to inactive the variants:
m1.variant(1).Active = 0;
m1.variant(2).Active = 0;
m1.variant(3).Active = 0;

% Then 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 = 300;
ncons = 35;
conc_vec = logspace(log10(0.2),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',[ 2.6 -1.6 0.83 0.0 ]);
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.028],...
    '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',[ 40 -2 0.83 0.01 ]);
% [Hill_2,GOF_2] = fit(X',Y',f,'StartPoint',[ 40 -2 0.83 0.01 ],...
%     'Lower',[30 -2.5 0.8 0.0],'Upper',[40 -1.5 0.84 0.02]);
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.028],...
    '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 2.7 mg, and an apparent Ki for the drug of 46
% micrograms/l. This is approaching the expected Ki = Kd = 1 micrograms/l.
% Notice the big difference of behavior with respect to the competitive
% inhibitor. Also, in this case we are getting Hill coefficients
% approaching 2, which is consistent with the inhibitor binding to both E
% and ES. We could get even better fits if we simulated two populations of
% inhibitor each with 1 binding site.

%% Developing a Drug Dosing Schedule
% Based on the ID50(ay 30 hours) of 2.7 mg we can try a schedule of 2
% mg/day. Counting on the progressive effect of the drug. Let's see what
% the stored repeated dosage looks like:

get(d2)

% First we set to inactive the variants:
m1.variant(1).Active = 0;
m1.variant(2).Active = 0;
m1.variant(3).Active = 0;

% Then 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');

d2.Amount = 2.0;
d2.Interval = 24;
d2.Rate = 2.0;
d2.RepeatCount = 25;

% We simulate for 960 hours (40 days). After 25 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 without 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 14]))
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 2 mg/24 hours for 25 days ');


%%
% Too much inhibition. We can probably achieve a very good result with a
% repeated smaller dose taken twice/day. 

d2.Amount = 0.68;
d2.Interval = 12;
d2.Rate = 0.68;
d2.RepeatCount = 50;

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 without 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 ))
legend('Aldosterone\_Plasma','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 0.68 mg/12 hours for 25 days ');


%%
% This looks great! We have achieved a stable normalization of the plasma
% aldosterone levels with only 2 daily doses and acceptably small daily
% fluctuations. However, it takes too long (200 hours) to achieve the
% normalization. We can combine a first large dose with smaller chasers to
% obtain a faster normalization.

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

d1.Amount = 2.0;
d2.Interval = 0;
d2.Rate = 2.0;
d2.RepeatCount = 0;
d2.StartTime = 6;

d2.Amount = 0.68;
d2.Interval = 12;
d2.Rate = 0.68;
d2.RepeatCount = 50;
d2.StartTime = 12;

% Initialize doses for analysis run.
doses = [d1, d2];

% Run simulation.
PKPD_Kinetics = sbiosimulate(m1, cs, [], doses);

%%
FIG_8 = figure;
    set(FIG_8,'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');
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))
legend('Aldosterone\_Plasma','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 single large dose + 0.68 mg/12 hours chase');

%% 
% Finally, we are interested in knowing whether the cc's for the
% aldosterone plasma concentration 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 500 hours (middle of a small fluctuation):

Stop = 500;
set(cs, 'StopTime', Stop);
delta = eps^(1/4);

nparams = size(params,1);
cc_Jk = zeros(1,nparams);
for n = 1:nparams 
    param_no = n; 
    k_ref = params(param_no).Value;
    k_range = linspace((1-delta)*k_ref,(1+delta)*k_ref,3);
    J_range = zeros(1,3);

    for i = 1:3
        params(param_no).Value = k_range(i);
        PKPD_sens = sbiosimulate(m1);
        ReactionFlux17 = (kcat_AS.Value*PKPD_sens.Data(end,11))*SG.capacity;
        J_range(i) = ReactionFlux17;
        C_range(i) = PKPD_sens.Data(end,4);                 
    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

%%
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

CONTROL_COEFF_2 = figure;
    set(CONTROL_COEFF_2,'Units','normalized','Position',[0.6 0.6 0.8 0.4],...
        'Name','Integrated Sensitivity Analysis');clf
    subplot(1,2,1)
        hbar = barh((cc_Jk(nz_ind)));
        haxes = get(hbar(1), 'Parent');
%         set(haxes,'XScale','log',...
%         'XMinorTick','on','XMinorGrid','on');
        set(haxes, 'ytick', 1:nnames);
        set(haxes, 'yticklabel', names);
        xlabel('cc(J\_ald,parameter)');
        ylabel('Parameters');
        set(gca,'XLim',[-1.2 1.2]);
        title('Flux Control Coefficients');

    subplot(1,2,2)
        hbar = barh((cc_Ck(nz_ind)));
        haxes = get(hbar(1), 'Parent');
%         set(haxes,'XScale','log',...
%         'XMinorTick','on','XMinorGrid','on');
        set(haxes, 'ytick', 1:nnames);
        set(haxes, 'yticklabel', names);
        xlabel('cc(C\_ald,parameter)');
        ylabel('Parameters');
        set(gca,'XLim',[-1.2 1.2]);
        title('Concentration Control Coefficients');
        
%% Conclusions
% Despite the Kd of the uncompetitive inhibitor for aldosterone synthase
% being 100 times higher than that of the competitive inhibitor, comparable
% levels of inhition are achieved with repeated doses of ~0.7 mg/12 hours
% instead of 30 mg/12 hours. Thus it can be expected that the side effects
% of the non-competitive drugs will be equally smaller.

%% General Rules:
% 
% # All competitive inhibitors have the same effects on substrate
% binding and catalysis. A competitive inhibitor will raise the apparent Km
% value for its substrate with no change in the apparent Vmax value. This
% characteristic will have physiological consequences on the observed
% efficacy of drugs. As an enzyme reaction is inhibited by a competitive
% inhibitor, there is an increase in the local concentration of substrate.
% As a consequence there is an increase in the apparent Ki value, and the
% inhibitor block is removed. Without a mechanism to clear the substrate, a
% competitive inhibitor will lose potency. This is not the case for a
% noncompetitive inhibitor.
% # A noncompetitive inhibitor binds equally well to both free enzyme
% and the enzyme-substrate complex. These binding events occur exclusively
% at a site distinct from the precise active site occupied by substrate. In
% contrast to a competitive inhibitor, a noncompetitive inhibitor will
% lower the apparent Vmax value, yet there is no effect on the apparent Km
% value for its substrate, and the Ki of the inhibitor does not change as a
% function of the substrate concentration. For this reason non-competitive
% inhibitor are usually far better drugs than competitive inhibitors.
% Unfortunately it is easier to find competitive than non- or uncompetitive
% drugs.
% # An uncompetitive inhibitor binds exclusively to the enzyme-substrate
% complex yielding an inactive enzyme-substrate-inhibitor complex: the
% apparent Vmax value and the apparent Km value should both decrease.
% Despite their rarity in drug discovery programs, uncompetitive inhibitors
% can be even better than non-competitive inhibitors. As the inhibitor
% decreases the enzyme activity, there is an increase in the local
% concentration of substrate. Without a mechanism to clear the buildup of
% substrate, the potency of the uncompetitive inhibitor actually increases.
% # The effectiveness of non-competitive (and uncompetitive) inhibitors
% is linked to the rate of degradation of the large pool of ESI complex
% that these inhibitors produce. The effect of the inhibitors is
% increased by high rates of degradation (short enzyme half-time), and
% diminished by low rates of degradation (long enzyme half-time).
%
%%
close all
