%% Exercise 5.4
% FRTF01 - Physiological Models and Computation, LTH, Lund University

close all
clear all

%-------------
% Parameters
%-------------
ke1 = 0.02;% min^-1
ke2 = 3e-4;% min^-1
ke3 = 8e-4;% min^-1
k12 = 0.1;% min^-1
k23 = 4e-3;% min^-1
k32 = 1e-3;% min^-1
VL = 2.7; %dl

%% Subproblem c

% Define system with oral input
%--------------
A = [-(ke1+k12) 0       0;...
    k12     -(ke2+k23)  k32;...
    0       k23     -(ke3+k32)];
B = [1 0 0]'; % Oral input
C = 1/VL*[0 0 1]; % Liver concentration [mg/dl]
D = [];
sys = ss(A,B,C,D);

% Oral dose 500 mg, assuming it takes 5 minutes to dissolve. 
% Make the trajectory over 168 hours
%----------------
u_oral = [100*ones(5,1); zeros(168*60-5,1)]; % 500 mg tablet dissolved over the first 5 minutes

% initial value
x0 = [0 0 0]'; %Assume we start with no drug in the body 

% simulation time
T = [0:1:(168*60-1)];

% Simulate
[y_1,T,x] = lsim(sys,u_oral,T,x0);

figure(1)
plot(T/60,y_1)
hold on
ylabel('Liver concentration [mg/dl]')
xlabel('Time [h]')

%% Subproblem d

% add the dosage in subproblem c every 24 hours 
u_onedose = [100*ones(5,1); zeros(24*60-5,1)];
u_oral24 = repmat(u_onedose,7,1); % 7 24-hour periods on 168 hours

% Simulate
[y_2,T,x] = lsim(sys,u_oral24,T,x0);

figure(1)
plot(T/60,y_2)

% alternative dosage: half the dosage every 12 hours
u_onedose = [100; 100; 50; zeros(12*60-3,1)];
u_oral12 = repmat(u_onedose,14,1); % 14 12-hour periods on 168 hours

% Simulate
[y_3,T,x] = lsim(sys,u_oral12,T,x0);

figure(1)
plot(T/60,y_3)


%% Subproblem e

% Define system with itravenous input 
%--------------
A = [-(ke1+k12) 0       0;...
    k12     -(ke2+k23)  k32;...
    0       k23     -(ke3+k32)];
B = [0 1 0]'; % IV input
C = 1/VL*[0 0 1]; % Liver concentration [mg/dl]
D = [];
sys = ss(A,B,C,D);

figure(2)
step(sys)

% From the step response it is evident that the static gain of the system
% is approx. 400 mg/dl. Could aso be calculated as follows: 
static_gain = dcgain(sys);

% If we want to achieve steady state liver concertration of 112 mg/dl we
% need the following magnitude of the constant iv-input: 
u_iv_mag = 112/static_gain; % U(0) = Y(0)/G(0);

% create the time trajectory of the IV-input
u_iv = u_iv_mag*ones(length(T),1);

% simulate
[y_4,T,x] = lsim(sys,u_iv,T,x0);

figure(1)
plot(T/60,y_4)
legend('500 mg once','500 mg every 24 h','250 mg every 12h','constant IV infusion',...
'Location','SouthEast')



