%%% Version: May 29th, 2017
%%%
%%% call by: Wajima_bloodCoagulation_PTtest_lowTF_MAIN
%%%
%%% Simulate the PT test (blood coagulation): tissue factor-fibrin 
%%% system and compute the corresponding input-response indices 
%%%
%%% Setting: Low initial tissue factor concentration of 5 pM
%%%
%%% by Wilhelm Huisinga and Jane Knoechel 
%%%

clear all; 

model.name = 'BloodCoagulationModel';

addpath(model.name)


%%% define tissue factor-fibrin system
%%%
I   = Wajima_bloodCoagulation_indexing;
I.neglectedStates = {}; 

X0  = Wajima_bloodCoagulation_initialValues(I);

%%% in the invitro setting Tmod was assumed to be zero Wajima (2009)
X0(I.Tmod) = 0;

%%% and the delution of the blood was considered by computing 1/3 of the
%%% orignal initial value
X0 = 1/3 * X0;

%%% assign variables to model structure
%%%

model.I        = I;
model.X0       = X0;

%%% the model can simulate different setting, where however the parameters
%%% are different (Wajima et al). In vivo all reactions are taking place,
%%% in vitro however the production of the factors is not active and thus
%%% assumed to be zero - this is realised by the string set in
%%% model.scenario
model.scenario = 'in_vitro'; % possible scenarios are in_vitro and in_vivo

%%% parameteres of the model were adapted by Gulati to simulate the influence of the
%%% snake venoms on the blood coagulation model. The original parameters
%%% are thus accessable by setting model.author to Wajima and the modified
%%% parameters by setting it to Gulati
model.author   = 'Wajima'; % possible authors are Wajima or Gulati

par = Wajima_bloodCoagulation_parameters(model);

%%% set parameter of intrinsic pathway to zero since otherwise during
%%% numerical solution of perturbed system with ode solver unphysiological
%%% activation of intrinsic pathway (set either v43 or v42 to zero)
%%%

par(I.v43) = 0;

%%% define input state
%%%

I.input  = I.TF;
u_ref = zeros(size(X0));
u_ref(I.input) = 5e-3;

%%% define output state
%%%
I.output = I.F;
model.h = @(x) x(:,I.output);

%%% assign variables to model structure
%%%
model.par   = par;
model.u_ref = u_ref;

%%% define system of ordinary differential equations
%%%
model.ode   = @(t,X) Wajima_bloodCoagulation_ode(t,X,par,model);

[~] = model_check(model);
%%% -----------------------------------------------------------------------
%%% solve the model for reference input
%%%

%%% the intesting dynamic for the PT test takes place one timespan of
%%% seconds 
SF_s_to_h = 1/3600;
SF_h_to_s = 1/SF_s_to_h;

model.SF_h_to_s = SF_h_to_s;

timespan = [0 240*SF_s_to_h];
X_init   = model.X0 + model.u_ref;

model.jac = @(t,X) Wajima_bloodCoagulation_jacobian(t,X,par,model);

options        = odeset('RelTol',1e-6,'AbsTol',1e-9,'NonNegative',1:I.nrOfStates,'Jacobian',model.jac);

[t_ref,x_ref]  = ode15s(model.ode,timespan,X_init,options);

%%% assign values to model structure
model.timespan = timespan;
model.X_init   = X_init;
model.t_ref    = t_ref;
model.x_ref    = x_ref;
model.solver.options = options;

%%% -----------------------------------------------------------------------
%%% graphical output for selected states
%%%
output.states = {'IXa','Xa','Fg','F','IIa','TF','VII_TF','VIIa_TF','Xa_Va','Va','VIII','VIIIa','IXa_VIIIa'};
output.color  = Wajima_bloodCoagulation_specifyColorStates(output.states);
output.legend = Wajima_bloodCoagulation_legendLabel(output.states);

model.output  = output;

FigNr = 1;
figure(FigNr); clf;
for i=1:length(output.states)
    plot(t_ref*SF_h_to_s,x_ref(:,I.(output.states{i})),'Color',output.color.(output.states{i}))
    hold on;
end
xlabel('t [s]'); ylabel('concentration [nM]');
xlim([-0 240]); ylim([1e-5 1e5])
axis square; set(gca,'yscale','log'); legend(output.legend,'Location','eastoutside');
fett(gcf); drawnow;

%%% =================================================================== %%%
%%% compute input-response indices
%%%
calculate_ir_indices = 1;

if calculate_ir_indices
       
    %%% perturb reference input by the following factors, i.e.,
    %%% u_per = u_ref * perturbationFactor
    model.perturbFac = 0.5:0.1:1.5;
        
    %%% check for consistency etc
    model.perturbFac(model.perturbFac==1) = []; % delete 0 entry, if it exists
    if sum(model.perturbFac<0)
        fprintf(' perturbation factors cannot be less than 0! ');
        error('')
    end;
    model.nrOfPerturb = length(model.perturbFac);
    
    %%% time span for which the input-response indices are computed
    %%% Note: the initial time is the start time of the trajectories
    %%% and the final time is the end time. The ir indices are by
    %%% definition zero for the end time, and except for the input, also 
    %%% zero for the start time
    tstarValues  = [0:(240*SF_s_to_h)/240:240*SF_s_to_h]; 
    model.tstarValues = tstarValues;

    start = tic;
    [ir]  = compute_input_response_indices_new(model);
    tictoc = toc(start); hours=floor(tictoc/(60*60)); minutes = floor(tictoc/60)-60*hours; seconds = round(tictoc-60*minutes-60*60*hours);
    fprintf('\n\n Total elapsed time: %d h %d min %2d sec',hours, minutes, seconds)

    
    identifier = 'lowTF';
    filename = sprintf('Wajima_bloodCoagulation_PTtest_%s',identifier);
    
    fprintf('\n\n results saved under filename = %s \n\n',filename);
    save(filename,'ir','model');
     
    %%% graphical output of the input-response indices for selected states
    %%% specified in output.states
    %%%
    FigNr = 2;
    linestyle = '-';
    if length(tstarValues)<5
           linestyle = '*';
    end;
    figure(FigNr); clf;
    for i=1:length(output.states)
        state = output.states{i};
        plot(tstarValues*SF_h_to_s,ir(:,I.(state)),linestyle,'Color',output.color.(state));
        hold on;
    end
    xlabel('t* [s]'); ylabel('ir(t*) index');
    xlim([0 240]); ylim([0 1e5]);
    axis square; set(gca,'yscale','lin'); legend(output.legend,'Location','eastoutside');
    fett(gcf)
    
end; % if calculate_ir_indices 


