%%% Version: September 11th, 2017
%%%
%%% call by: Create_figures_Knoechel_etal_InputReponseIndex
%%%
%%% with this Main function you can either recreate a specific or all
%%% figures based on model simulations given in the following article:
%%%
%%% Note in order to recreate figures as in publication do not change the
%%% switch environment
%%%
%%% Knoechel, Kloft and Huisinga, "Reducing complex
%%% systems pharmacology models based on a novel input-response Index"
%%% J Pharmacokinet Pharmacodyn, 2017 (submitted)
%%%
%%% by Jane Knoechel
%%%

%%% specify the figure names to be recreated (3a and 3b,5,6,7,9a and 9b)
%%%
figName = {'3a' '3b' '5' '6' '7' '9a' '9b' '12' '13' '14'};

%%% create chosen figure by calling necessary matlab functions
%%%
model.name = 'BloodCoagulationModel';

addpath(model.name)

if find(ismember('3a',figName),1,'first')
    %%% load the input-response indices and variables for model specification
    %%%
    file    = dir('Wajima_bloodCoagulation_snakeVenom_1h.mat');
    load(file.name)
    
    output.states={'P','Fg','IIa','APC_PS','IIa_Tmod','APC','AVenom','CVenom'};
    output.color  = Wajima_bloodCoagulation_specifyColorStates(output.states);
    output.legend = Wajima_bloodCoagulation_legendLabel(output.states);
    %%% figure
    %%%
    figure(31); clf;
    for i=1:length(output.states)
        state = output.states{i};
        plot(model.tstarValues,ir(:,model.I.(state)),'Color',output.color.(state));
        hold on;
    end
    xlabel('Time* [h]'); ylabel('Input-response index');
    xlim([0 0.15]); ylim([0 7e4]);
    axis square; set(gca,'yscale','lin'); legend(output.legend,'Location','eastoutside');
    fett(gcf)
    %%% inset
    %%%
    axes('Position',[0.461360123647604 0.583314794215795 0.2 0.2])
    box on
    for i=1:length(output.states)
        state = output.states{i};
        plot(model.tstarValues,ir(:,model.I.(state)),'Color',output.color.(state))
        hold on;
    end
    set(gca,'XLim',[0 model.timespan(end)])
    set(gca,'YLim',[0 7e4])
    axis square
    fett(gcf)

end
if find(ismember('3b',figName),1,'first')
    %%% load the input-response indices and variables for model specification
    %%%
    file    = dir('Wajima_bloodCoagulation_snakeVenom_1h.mat');
    load(file.name)
    
    output.states={'P','Fg','IIa','APC_PS','IIa_Tmod','APC','AVenom','CVenom','XIa','IXa','IXa_VIIIa','VIIIa','Xa','V','Va','Xa_Va','II','Tmod'};
    output.color  = Wajima_bloodCoagulation_specifyColorStates(output.states);
    output.legend = Wajima_bloodCoagulation_legendLabel(output.states);
   
    %%% figure
    %%%
    figure(32); clf;
    for i=1:length(output.states)
        state = output.states{i};
        plot(model.t_ref,model.x_ref(:,model.I.(state)),'Color',output.color.(state));
        hold on;
    end
    xlabel('Time [h]'); ylabel('Concentration [nM]');
    xlim([0 0.15]); ylim([1e-5 1e6]);
    axis square; set(gca,'yscale','log'); legend(output.legend,'Location','eastoutside');
    fett(gcf)
    %%% inset
    %%%
    axes('Position',[0.461360123647604 0.583314794215795 0.2 0.2])
    box on
    for i=1:length(output.states)
        state = output.states{i};
        plot(model.t_ref,model.x_ref(:,model.I.(state)),'Color',output.color.(state))
        hold on;
    end
    set(gca,'XLim',[0 model.timespan(end)])
    set(gca,'YLim',[1e-5 1e6])
    axis square;set(gca,'yscale','log');
    fett(gcf)
end

if find(ismember('6',figName),1,'first')
    %%% load the input-response indices and variables for model specification
    %%%
    file    = dir('Wajima_bloodCoagulation_snakeVenom_1h.mat');
    load(file.name)
    
    I = model.I;
    %%% compute rates of degradation for fibrinogen
    %%%
    v_P = (model.par(I.v15)*model.x_ref(:,I.P).*model.x_ref(:,I.Fg)./(model.par(I.k15)+model.x_ref(:,I.P)));
    v_IIa = model.par(I.v14)*model.x_ref(:,I.IIa).*model.x_ref(:,I.Fg)./(model.par(I.k14)+model.x_ref(:,I.IIa));
    v_kdegFg =model.par(I.degFg)*model.x_ref(:,I.Fg);
    
    figure1=figure(6);
    plot(model.t_ref,v_P,'k',model.t_ref,v_IIa,'-.k',model.t_ref,v_kdegFg,'--k')
    hold on;
    % Create textarrow
    annotation(figure1,'textarrow',[0.375563872614212 0.343879174005247],...
        [0.752678013622442 0.692055099273165],'String',{'v_{P}'},'LineWidth',3,...
        'FontWeight','bold',...
        'FontSize',20);
    
    % Create textarrow
    annotation(figure1,'textarrow',[0.420428812767334 0.37375184213364],...
        [0.300112758468313 0.169300745120148],'String',{'v_{IIa}'},'LineWidth',3,...
        'FontWeight','bold',...
        'FontSize',20);
    
    % Create textarrow
    annotation(figure1,'textarrow',[0.29796511627907 0.294651522231408],...
        [0.167808219178082 0.109522376460908],'String',{'v_{kdeg Fg}'},...
        'LineWidth',3,...
        'FontWeight','bold',...
        'FontSize',20);
    
    xlabel('Time [h]','FontSize',20);
    ylabel('Transformation rate of fibrinogen [nM/h]','FontSize',20);
    legend({'v_{P}','v_{IIa}','v_{kdeg Fg}'})
    legend('hide')
    xlim([0 0.3])
    axis square
    fett(gcf)
    %%% inset
    %%%
    axes('Parent',figure1,'Position',[0.600604309694116 0.645643561339083 0.2 0.2])
    box on
    plot(model.t_ref,v_P,'k',model.t_ref,v_IIa,'-.k',model.t_ref,v_kdegFg*ones(1,length(model.t_ref)),'--k')
    xlim([0 1])
    axis square
    fett(gcf)
end
if find(ismember('7',figName),1,'first')
    %%% load the input-response indices and variables for model specification
    %%%
    file    = dir('Wajima_bloodCoagulation_snakeVenom_40h.mat');
    load(file.name)
    
    output = model.output;
    %%% figure
    %%%
    figure(7); clf;
    for i=1:length(output.states)
        state = output.states{i};
        plot(model.tstarValues,ir(:,I.(state)),'Color',output.color.(state));
        hold on;
    end
    xlabel('Time* [h]'); ylabel('Input-response index');
    xlim([0 0.15]); ylim([0 1.5e3]);
    axis square; set(gca,'yscale','lin'); legend(output.legend,'Location','eastoutside');
    fett(gcf)
    %%% inset
    %%%
    axes('Position',[0.461360123647604 0.583314794215795 0.2 0.2])
    box on
    for i=1:length(output.states)
        state = output.states{i};
        plot(model.tstarValues,ir(:,I.(state)),'Color',output.color.(state))
        hold on;
    end
    set(gca,'XLim',[0 model.timespan(end)])
    set(gca,'YLim',[0 1.5e3])
    axis square
    fett(gcf)
end
if find(ismember('9a',figName),1,'first')
    %%% load the input-response indices and variables for model specification
    %%%
    file    = dir('Wajima_bloodCoagulation_PTtest_lowTF.mat');
    load(file.name)
    
    output = model.output;
    I      = model.I; 
    %%% figure
    %%%
    figure(91); clf;
    for i=1:length(output.states)
        state = output.states{i};
        plot(model.tstarValues*model.SF_h_to_s,ir(:,I.(state)),'-','Color',output.color.(state));
        hold on;
    end
    xlabel('Time* [s]'); ylabel('Input-response index');
    xlim([0 240]); ylim([0 12e5]);
    axis square; set(gca,'yscale','lin'); legend(output.legend,'Location','eastoutside');
    fett(gcf)
end
if find(ismember('9b',figName),1,'first')
    %%% load the input-response indices and variables for model specification
    %%%
    file    = dir('Wajima_bloodCoagulation_PTtest_highTF.mat');
    load(file.name)
    
    output = model.output;
    %%% figure
    %%%
    figure(92); clf;
    for i=1:length(output.states)
        state = output.states{i};
        plot(model.tstarValues*model.SF_h_to_s,ir(:,I.(state)),'-','Color',output.color.(state));
        hold on;
    end
    xlabel('Time* [s]'); ylabel('Input-response index');
    xlim([0 30]); ylim([0 70]);
    axis square; set(gca,'yscale','lin'); legend(output.legend,'Location','eastoutside');
    fett(gcf)
        
end
if find(ismember('12',figName),1,'first')
    file = dir('Wajima_bloodCoagulation_snakeVenom_40h.mat');
    load(file.name)
    %%% output for 62 state variable model
    t_ref = model.t_ref;   
    x_ref = model.x_ref;
    y_ref = model.h(x_ref);
    I     = model.I;
    
    y_refcolor  = Wajima_bloodCoagulation_specifyColorStates({'Fg'});
    
    par = model.par;
    %%% compute knock-outs 
    
    odeFcn   = @(t,X,par) Wajima_bloodCoagulation_ode(t,X,par,model);
    %%% knockout of X and V
    par_knockout_X_V = par;
    par_knockout_X_V(I.aX)=0;
    par_knockout_X_V(I.pV)=0;

    X_init_knockout_X_V = model.X_init;
    X_init_knockout_X_V(I.X) = 0;
    X_init_knockout_X_V(I.V) = 0;
    
    [~,x_knockout_X_V]  = ode15s(@(t,X) odeFcn(t,X,par_knockout_X_V),t_ref,X_init_knockout_X_V,model.solver.options);
    y_knockout_X_V = model.h(x_knockout_X_V);
    
    %%% knockout of V
    par_knockout_V = par;
    par_knockout_V(I.pV)=0;

    X_init_knockout_V = model.X_init;
    X_init_knockout_V(I.V) = 0;
    
    [~,x_knockout_V]  = ode15s(@(t,X) odeFcn(t,X,par_knockout_V),t_ref,X_init_knockout_V,model.solver.options);
    y_knockout_V = model.h(x_knockout_V);
    
    %%% knockout of Pg
    par_knockout_Pg = par;
    par_knockout_Pg(I.pPg)=0;
    
    X_init_knockout_Pg = model.X_init;
    X_init_knockout_Pg(I.Pg) = 0;
    
    [~,x_knockout_Pg]  = ode15s(@(t,X) odeFcn(t,X,par_knockout_Pg),t_ref,X_init_knockout_Pg,model.solver.options);
    y_knockout_Pg = model.h(x_knockout_Pg);
    
    SF_nm_to_um =1/1000;
    Outputlegend = {'62-state variable model','knock-out X & V','knock-out V','knock-out Pg'};
    
    figure(12); clf;
    plot(t_ref,y_ref*SF_nm_to_um,'Color',y_refcolor.Fg);
    hold on;
    plot(t_ref,y_knockout_X_V*SF_nm_to_um,'-.','Color',y_refcolor.Fg);
    plot(t_ref,y_knockout_V*SF_nm_to_um,'--','Color',y_refcolor.Fg);
    plot(t_ref,y_knockout_Pg*SF_nm_to_um,':','Color',y_refcolor.Fg);
    
    xlabel('Time [h]'); ylabel('Fibrinogen concentration [uM]');
    xlim([0 10]); ylim([1e-3 1e1]);
    axis square; set(gca,'yscale','log'); legend(Outputlegend,'Location','eastoutside');
    fett(gcf)
    %%% inset
    %%%
    axes('Position',[0.461360123647604 0.583314794215795 0.2 0.2])
    box on
    plot(t_ref,y_ref*SF_nm_to_um,'Color',y_refcolor.Fg);
    hold on;
    plot(t_ref,y_knockout_X_V*SF_nm_to_um,'-.','Color',y_refcolor.Fg);
    plot(t_ref,y_knockout_V*SF_nm_to_um,'--','Color',y_refcolor.Fg);
    plot(t_ref,y_knockout_Pg*SF_nm_to_um,':','Color',y_refcolor.Fg);
    
    set(gca,'XLim',[0 1])
    axis square;set(gca,'yscale','lin');
    fett(gcf)
end
if find(ismember('13',figName),1,'first')
   %%% load the input-response indices and variables for model specification
    %%%
    file    = dir('Wajima_bloodCoagulation_snakeVenom_1h.mat');
    load(file.name)
    
    output.states = setdiff(model.I.stateName,{'P','Fg','IIa','APC_PS','IIa_Tmod','APC','AVenom','CVenom'});
    output.color  = Wajima_bloodCoagulation_specifyColorStates(output.states);
    output.legend = Wajima_bloodCoagulation_legendLabel(output.states);
   
    %%% figure
    %%%
    figure(13); clf;
    for i=1:length(output.states)
        state = output.states{i};
        plot(model.t_ref,model.x_ref(:,model.I.(state)),'Color',output.color.(state));
        hold on;
    end
    xlabel('Time [h]'); ylabel('Concentration [nM]');
    xlim([0 0.15]); ylim([1e-5 1e6]);
    axis square; set(gca,'yscale','log'); legend(output.legend,'Location','eastoutside');
    fett(gcf)
    %%% inset
    %%%
    axes('Position',[0.461360123647604 0.583314794215795 0.2 0.2])
    box on
    for i=1:length(output.states)
        state = output.states{i};
        plot(model.t_ref,model.x_ref(:,model.I.(state)),'Color',output.color.(state))
        hold on;
    end
    set(gca,'XLim',[0 model.timespan(end)])
    set(gca,'YLim',[1e-5 1e6])
    axis square;set(gca,'yscale','log');
    fett(gcf)
    %%% compute legend for figure
    figure(131); clf;
    ax1 = subplot(1,2,1);
    for i=1:length(output.states)/2
        state = output.states{i};
        p=plot(model.t_ref,model.x_ref(:,model.I.(state)),'Color',output.color.(state));
        hold on;
        set(p,'visible','off')
    end
    xlabel('Time [h]'); ylabel('Concentration [nM]');
    xlim([0 0.15]); ylim([1e-5 1e6]);
    set(gca,'visible','off')
    axis square; set(gca,'yscale','log'); legend(output.legend(1:length(output.states)/2),'Location','eastoutside');
    fett(gcf)
    ax1 = subplot(1,2,2);
    for i=length(output.states)/2+1:length(output.states)
        state = output.states{i};
        p=plot(model.t_ref,model.x_ref(:,model.I.(state)),'Color',output.color.(state));
        hold on;
        set(p,'visible','off')
    end
    xlabel('Time [h]'); ylabel('Concentration [nM]');
    xlim([0 0.15]); ylim([1e-5 1e6]);
    set(gca,'visible','off')
    axis square; set(gca,'yscale','log'); legend(output.legend(length(output.states)/2+1:length(output.states)),'Location','eastoutside');
    fett(gcf)
end
if find(ismember('14',figName),1,'first')
    %%% compute comparison of all state variables except output between 62
    %%% state variable model and elimination-reduced model
    
    
    %%% load the input-response indices and variables for model specification
    %%%
    file    = dir('Wajima_bloodCoagulation_snakeVenom_1h.mat');
    load(file.name)
    
    %%% output for 62 state variable model
    t_ref   = model.t_ref;   
    x_ref   = model.x_ref;
    I       = model.I;
    X_init  = model.X_init;
    
    redmodel = model;
    
    odeFcn               = @(t,X,model)  Wajima_bloodCoagulation_ode(t,X,model.par,model);
    
    redmodel.I.environStates   = {'Pg','II','PC','PS','Tmod'};
    redmodel.I.selectedStates  = {'Pg','P','Fg','II','IIa','APC_PS','Tmod','IIa_Tmod','PC','APC','PS','AVenom','CVenom'};
    dynStates = {'P','IIa','APC_PS','IIa_Tmod','APC','AVenom','CVenom'};
    
    %%% create transformation matrix
    %%%
    redmodel.transf = create_transf_matrices(redmodel,'Galerkin');
    
    %%% set initial value for reduced model
    %%%
    X_init_red = redmodel.transf.M*X_init;
    
    %%% reset solver options
    redmodel.solver.options.NonNegative = 1:length(redmodel.I.selectedStates);
    redmodel.solver.options.Jacobian     = [];
    
    %%% compute solution of reduced model with neglected states
    %%%
    [~,x_tmp] = ode15s(@(t,X) odeFcn(t,X,redmodel),t_ref,X_init_red,redmodel.solver.options);
    
    %%%  map reduced solution back to full original system to compute the output
    %%%
    x_red      = x_tmp*redmodel.transf.invM';
    
    dynStatesColor  = Wajima_bloodCoagulation_specifyColorStates(dynStates);
    dynStatesLegend = Wajima_bloodCoagulation_legendLabel(dynStates);
    
    for i = 1:length(dynStates)
        LegendDynStatesFullAndRed(i*2-1) = dynStatesLegend(i);
        LegendDynStatesFullAndRed(i*2)   = strcat(dynStatesLegend(i),'_{red}');
    end
    
    figure(14)
    for i=1:length(dynStates)
        plot(t_ref,x_ref(:,model.I.(dynStates{i})),'color',dynStatesColor.(dynStates{i}))
        hold on;
        plot(t_ref,x_red(:,model.I.(dynStates{i})),'--','color',dynStatesColor.(dynStates{i}))
    end
    xlabel('Time [h]'); ylabel('Concentration [nM]');
    xlim([-0 0.15]); ylim([1e-5 1e6])
    axis square; set(gca,'yscale','log'); legend(LegendDynStatesFullAndRed,'Location','eastoutside');
    fett(gcf); drawnow;
    
    axes('Position',[0.461360123647604 0.583314794215795 0.2 0.2])
    box on
    for i=1:length(dynStates)
        plot(t_ref,x_ref(:,model.I.(dynStates{i})),'color',dynStatesColor.(dynStates{i}))
        hold on;
        plot(t_ref,x_red(:,model.I.(dynStates{i})),'--','color',dynStatesColor.(dynStates{i}))
    end
    set(gca,'yscale','log');
    set(gca,'XLim',[-0 1])
    set(gca,'YLim',[1e-5 1e6])
    axis square
    fett(gcf)
end

if find(ismember('5',figName),1,'first')
    %%% figure
    %%%
    figure(5); clf;
    run Wajima_bloodCoagulation_snakeVenom_MOR_40h
end