% 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 (non-competitive drug)

% This example shows how to analyze a pharmacokinetics-pharmacodynamics
% model in which a drug, acting as a non-competitive 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 a non-competitive 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('PKPDnc');
% addpath(genpath('vhline'));

% or alternatively:

% m1 = sbmlimport('PKPDnc.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 = 'PKPDnc.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
% non-competitive inhibitor of this enzyme.

%%
% <<../PKPDnc.jpg>>

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

%%
% <<../Non-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:
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_E = sbioselect(m1,'Name','kon_E');
koff_E = sbioselect(m1,'Name','koff_E');
kcat_E = sbioselect(m1,'Name','kcat_E');
kon_EI = sbioselect(m1,'Name','kon_EI');
koff_EI = sbioselect(m1,'Name','koff_EI');
kon_ES = sbioselect(m1,'Name','kon_ES');
koff_ES = sbioselect(m1,'Name','koff_ES');
kon_ESI = sbioselect(m1,'Name','kon_ESI');
koff_ESI = sbioselect(m1,'Name','koff_ESI');
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
kdeg_ESI = sbioselect(m1,'Name','kdeg_ESI');    % degradation ESI complex
kti_ald = sbioselect(m1,'Name','kti_ald');
kit_ald = sbioselect(m1,'Name','kit_ald');

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

Drug_Kd = koff_EI.Value/kon_EI.Value;
display(['Drug_Kd = ' num2str(Drug_Kd) ' microgram/L']);

% The Kd for this non-competitive 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

% Please notice that states are logged in the following order (DataNames);

%    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 
% To make sure nothing has changed in the model in the absence of drug we
% repeat the sensitivity analysis.

% Configure task specific stop time. We simulate up to 30 hours.

set(cs, 'StopTime', 30.0);  
params = sbioselect(m1,'Type','parameter');

% 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);
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);
        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 = diff(J_range/J_range(2))./diff(k_range/k_range(2)); 
%     dCdk = diff(C_range/C_range(2))./diff(k_range/k_range(2));
    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) = mean(dJdk);  
    cc_Ck(n) = mean(dCdk);  
    params(param_no).Value = k_ref;
end

%%
% and derive the value of each cc

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.

%% Dose Response Curves
% Let's also check that the amount of aldosterone synthase that needs to be
% inhited has not changed because of the additional equilibria included in
% this model.

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.1],...
    'String',{[string1 string2 string3]},...
    'FontWeight','bold',...
    'FontSize',14,...
    'FitBoxToText','off',...
    'BackgroundColor',[1 0.968627452850342 0.921568632125854],...
    'Color',[1 0 0]);

%%
% Everything look OK, and we have confirmed that we need to bring the activity of
% aldosterone synthase to ~ 5% 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;

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

%--------------------------------------------------------------------------
% We use an 8 parameter logistic fit:
         
f = fittype('(c - d)/(1 + (a/x)^b) + (g - h)/(1 + (e/x)^f)');
[Hill_1B,GOF_1B] = fit(X',Y',f,'StartPoint',[ 1.9 -1.8 0.83 0.4 4.2 -1.2 0.4 0.0 ]);
Hill_1B
%--------------------------------------------------------------------------


semilogx(Dose,Hill_1(Dose)','-r',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','y',...
             'MarkerSize',5,'Parent',axes_dr1);

semilogx(Dose,Hill_1B(Dose)','-g',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','y',...
             'MarkerSize',5,'Parent',axes_dr1);
         
legend('Plasma Ald','4 params Logistic Fit','8 params Logistic Fit',...
    'Location','NorthEast');
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',[ 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

%--------------------------------------------------------------------------
% We use an 8 parameter logistic fit:
f = fittype('(c - d)/(1 + (a/x)^b) + (g - h)/(1 + (e/x)^f)');
[Hill_2B,GOF_2B] = fit(X',Y',f,'StartPoint',[ 30 -2 0.83 0.4 80 -1 0.4 0]);
Hill_2B
%--------------------------------------------------------------------------

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

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

semilogx(Drug_conc,Hill_2B(Drug_conc)','-m',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','y',...
             'MarkerSize',5,'Parent',axes_dr2);
         
legend('Plasma Ald','4 params Logistic Fit','8 params Logistic Fit',...
    'Location','NorthEast');
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 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 can get an even better fit if we simulate two populations of
% binding proteins (ES and E) each with 'n' binding site. This fit is shown
% with the green (top figure) and magenta lines (bottom figure). However,
% as far as the ID50 and the Ki are concerned the simpler fit is
% sufficient.

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

get(d2)

% It look reasonable.
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 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 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.62;
d2.Interval = 12;
d2.Rate = 0.62;
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 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 0.62 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;
d1.Interval = 0;
d1.Rate = 2.0;
d1.RepeatCount = 0;
d1.StartTime = 6;

d2.Amount = 0.62;
d2.Interval = 12;
d2.Rate = 0.62;
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 6 15]))
legend('Aldosterone\_Plasma','Aldosterone\_IF','Aldosterone\_OT',...
    'Location','Best')
hline([0.34 0.34],{'g-','LineWidth',2});
% 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.62 mg/12 hours chase');

%% 
% 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 500 hours
% (middle of a small fluctuation):

Stop = 500;
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, [], doses);
        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 = diff(J_range/J_range(2))./diff(k_range/k_range(2)); 
%     dCdk = diff(C_range/C_range(2))./diff(k_range/k_range(2)); 
    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) = mean(dJdk);  
    cc_Ck(n) = mean(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.4));
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.4 0.4]);

% 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.4));
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.4 0.4]);
        
%%
% The presence of the inhibitor affects the cc's for the synthesis of
% aldosterone in a completely different way from the competitive inhibitor.
% One interesting thing to notice is that, besides the obvious kon_ES and
% koff_ES of the inhibitor for the ES complex, and the parameters that
% affect the distribution of drug and aldosterone between compartments, the
% two parameters that produce the most significant inhibition of
% aldosterone synthesis are the rate of degradation of the ESI and ES
% complexes (kdeg_ES and kdeg_ESI). This is because, as new enzyme is being
% synthesized, the pool of the ESI complex becomes progressively larger,
% and at some point it starts increasing the ES pool producing an increase
% of the rate of aldosterone production. On the other hand the two
% parameters whose increase produces the largest increase in aldosterone
% synthesis and plasma aldosterone concentration are kcat_E and kdeg_E.
% Thus, the half-life of the enzyme affects the plasma aldosterone level in
% a complex way. We will now investigate this issue. We recall here that
% the half-life can be calculated from the degradation rate constants as:
% t_half = ln(2)/kdeg.

%% An allosteric drug that increases the rate of degradation of EI and ESI.
% Based on the previous analysis, we expect that if a modified drug
% increases the rate of degradation of the ESI complex by virtue of an
% allosteric effect, plasma aldosterone levels will drop further. Let's
% check that by increasing by 10 the relevant rate constants kdeg_EI and 
% kdeg_ESI.

kdeg_EI_ref = kdeg_EI.Value;
kdeg_ESI_ref = kdeg_ESI.Value;
kdeg_E_ref = kdeg_E.Value;
kdeg_ES_ref = kdeg_ES.Value;

kdeg_EI.Value = 2;
kdeg_ESI.Value = 2;
kdeg_E.Value = 0.2;
kdeg_ES.Value = 0.2;

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

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

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

hold(axesDal,'all');
box(axesDal,'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 allosteric drug that shortens 10 times EI and ESI half-life');
hold on

%%
% Yes, our expectation is confirmed. However, if the half-life of all the
% enzyme forms (also E and ES) is shortened the effect is completely different:

kdeg_E.Value = 2;
kdeg_ES.Value = 2;

% Run simulation.
PKPD_Kinetics = sbiosimulate(m1, cs, [], doses);
plot(PKPD_Kinetics.Time,PKPD_Kinetics.Data(:,[4]),'-c','Linewidth',1)
legend('Aldosterone\_Plasma','Aldosterone\_IF',...
    'Aldosterone\_OT','Aldosterone\_Plasma 2','Location','Best')

%% Variants with altered half-life of aldosterone synthase.
% Thus, the dominant effect is associated with the half-life of the
% free-enzyme (kdeg_E). This situation can be equated to that of a mutant
% clone of the surrenal tumor cells that has developed a high rate of
% synthesis and degradation of aldosterone synthase. The result is
% that the half-life of the enzyme is significantly shortened, although the
% steady state amount of the enzyme in the cell remains unchanged. Under
% these conditions, the dosing schedule we had developed is completely
% ineffective. Is there any way we can improve our drug in order to restore
% its efficacy? For this purpose, let's see if the control coefficients of
% the system parameters  with respect to plasma aldosterone have changed:

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

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

for n = 1:nparams 
    param_no = n; 
    k_ref = params(param_no).Value;
    delta = 0.009;
    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, [], doses);
        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 = diff(J_range/J_range(2))./diff(k_range/k_range(2)); 
    dCdk = diff(C_range/C_range(2))./diff(k_range/k_range(2)); 
    cc_Jk(n) = mean(dJdk);  
    cc_Ck(n) = mean(dCdk);  
    params(param_no).Value = k_ref;
end

%%
% and derive the value of each cc

CONTROL_COEFF_3 = figure;
    set(CONTROL_COEFF_3,'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_3,...
    '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.02));
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_3,...
    '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.02 0.02]);

% 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_3,...
    '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.02));
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_3,...
    '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.02 0.02]);
        
%%
% We can see that now the rate of degradation of EI (kdeg_EI) or ESI
% (kdeg_ESI) is no longer an important parameter. The best way we could
% modify our drug is by increasing the kon for the inhibitor
% (kon_ES,kon_EI). This is the same as decreasing the Kd. Let's see what
% happens if we do that:

kon_EI_ref = kon_EI.Value;
koff_EI_ref = koff_EI.Value;
kon_ES_ref = kon_ES.Value;
koff_ES_ref = koff_ES.Value;

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

kon_EI.Value = 0.65;
kon_ES.Value = 0.65;

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

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

hold(axesD2,'all');
box(axesD2,'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 high affinity non-competitive drug ');

%%
% We would have to increase the kon for the drug by 6-7 times in order to
% restore its effect. The alternative is to increase the drug dosing.
% First, we restore the original values for kon and koff for ES and EI.

koff_EI.Value = 0.1;
kon_EI.Value = 0.1;
koff_ES.Value = 0.1;
kon_ES.Value = 0.1;

%%
% Then we set the new dosing:

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


d1.Amount = 4.2;
d1.Interval = 0;
d1.Rate = 4.2;
d1.RepeatCount = 0;
d1.StartTime = 6;

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

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

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

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

hold(axesD3,'all');
box(axesD3,'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 higher dose drug');

% Thus, a decrease in the half-life of the enzyme can be corrected by
% either a comparable decrease in the Kd of the inhibitor for the enzyme,
% or by a comparable increase in the inhibitor dose.

%% Conclusions
% Although the Kd of the non-competitive inhibitor for aldosterone synthase
% is 100 times higher than that of the competitive inhibitor, comparable
% levels of inhition are achieved with repeated doses of ~0.6 mg/12 hours
% instead of 25 mg/12 hours. Thus, it is possible that the side effects
% of the non-competitive drugs will be comparably 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. We have observed this increase in the
% apparent Ki value when we determined a dose-response curve for the plasma
% aldosterone level with respect to the concentration of drug in the
% surrenal gland. 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. In fact we have observed that
% the apparent Ki value determined from the dose-response curve for the
% plasma aldosterone level with respect to the concentration of drug in the
% surrenal gland approached the real Kd value of the drug for either E or
% ES. 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 half-life of the enzyme. Smaller half-lifes require either
% a higher dose of the inhibitor, or an inhibitor with a comparably smaller
% Kd for the enzyme.
%
%%
close all
% save(savefile);
