% 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.
%
%% Transport of a solute between two compartments

% This example explains the application of the mass action law to the  
% transfer of a solute between two compartments.

%%
% We start by loading the project :
clear, clc

sbioloadproject('../TOOLBOXES/PKPD/Two_boxes');

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

V1 = sbioselect(m1,'Name','V1');
V2 = sbioselect(m1,'Name','V2');
A1 = sbioselect(m1,'Name','A1');
A2 = sbioselect(m1,'Name','A2');
kV1V2 = sbioselect(m1,'Name','kV1V2');
kV2V1 = sbioselect(m1,'Name','kV2V1');

% We also store some initial values:
A1_Init = A1.InitialAmount;
A2_Init = A2.InitialAmount;

% And we set the rate constants
kV1V2.Value = 0.01;
kV2V1.Value = 0.01;

%% 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:
Stop = 300;
set(cs, 'SolverType', 'ode15s');
set(cs, 'StopTime', Stop);
set(cs, 'TimeUnit','second');
set(cs.RunTimeOptions, 'StatesToLog', 'All');
set(cs.CompileOptions, 'UnitConversion', true);
set(cs.CompileOptions, 'DimensionalAnalysis', true);
cs
%%

Two_boxes_kinetics = sbiosimulate(m1, cs);
Two_boxes_kinetics.DataNames
m1.species


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

hold(axes1,'all');
box(axes1,'on');
plot(Two_boxes_kinetics.Time,Two_boxes_kinetics.Data(:,:))
legend('A1','A2')
% set(gca,'YLim',[-0.05 1.05]);
% set(gca,'XLim',[-0.00001 1.05*Stop]);
ylabel('Species Concentration (moles/liter)')
xlabel('Time (seconds)')
title('Equilibrium between two boxes ');

%% Scaling of rate constants to achieve equilibrium

A1.InitialAmount = A1_Init ;
A2.InitialAmount = A2_Init;
kV1V2.Value = 1.0;
kV2V1.Value = 0.01;

Two_boxes_kinetics = sbiosimulate(m1, cs);

FIG_2 = figure;set(FIG_2,'Units','normalized','Position',[0.6 0.6 0.4 0.8]);clf
subplot(2,1,1)
plot(Two_boxes_kinetics.Time,Two_boxes_kinetics.Data(:,1),'-b')
legend('A1');ylabel('Species Concentration (moles/liter)');xlabel('Time (seconds)')
title('Equilibrium between two boxes ');
subplot(2,1,2)
plot(Two_boxes_kinetics.Time,Two_boxes_kinetics.Data(:,2),'-r')
legend('A2');ylabel('Species Concentration (moles/liter)');xlabel('Time (seconds)')
