function[biomarks] = calculate10biomark(s,model)
% input
% s: struct with simulation results
% model: string indicating 'ORdmD' or 'TorORd'
%
% matrix with nx10 biomarkers, being n the numbers of models simulated in
% s.
% biomarks(j,:)=[apd90 tri90_30 tri90_50 qnet cai_syst cai_diast ctd90
% ctd50 EMw EAD];
    

biomarks=zeros(size(s.result,1),10); 
for j=1:size(s.result,1) 
    for i=1:size(s.result{j,1},1)
        t=s.result{j,1}{i,1};
        v=s.result{j,2}{i,1}(:,1);
        cai=s.result{j,2}{i,1}(:,6)*10^3;
        [apd90(i),~,dvdtmax(i),rmp(i),vpeak(i), EAD(i)]=APD(t,v,90);
        apd50(i)=APD(t,v,50);
        apd30(i)=APD(t,v,30);
        tri90_50(i)=apd90(i)-apd50(i);
        tri90_30(i)=apd90(i)-apd30(i);

        [ctd90(i),cai_syst(i),cai_diast(i),EAD2(i)]=CTD(t,cai,90, EAD(i));
        ctd50(i)=CTD(t,cai,50, EAD(i));

        EMw(i)=ctd90(i)-apd90(i);

        if strcmp(model,'ORdmJ')
            [qnet(i), qna(i), qnal(i), qcal(i), qto(i), qkr(i), qks(i), qk1(i)] = calculateCurrentsIntegrals(s.result{j,1},s.result{j,3},model,i);
        elseif strcmp(model,'TorORd')
            s.params(j).model=@model_Torord_dynCl;
            currents = getCurrentsStructure(s.result{j,1}, s.result{j,2}, s.params(j), 0);
            [qnet(i), qna(i), qnal(i), qcal(i), qto(i), qkr(i), qks(i), qk1(i)] = calculateCurrentsIntegrals(s.result{j,1},currents,model,i);
        end
    end
    if sum(EAD ==0)==3 && sum(EAD2 ==0)==3
        biomarks(j,:)=[mean(apd90,2) mean(tri90_30,2) mean(tri90_50,2) mean(qnet,2) ...
            mean(cai_syst,2) mean(cai_diast,2) mean(ctd90,2) mean(ctd50,2) ...
            mean(EMw,2) mean(EAD,2)];
    else
        biomarks(j,1:17)=NaN;
    end
    clearvars EAD EAD2
end
end

function [APDxx,Vpp,dVdt_max,Vr,Vp, EAD]=APD(t,v,xx)
% APD at different repolarization % (xx, eg 90%)
% t: vector time
% v: vector potential (1 AP)
% xx: repolarization %
% Other parameters: amplitude (Vpp), max upstroke (dVdt_max), Vrest, Vpeak
% Ex.: apd90=APD(t,v,90);  

[Vp, Vp_index]=max(v); %Caution with elevated domes
[Vr, Vr_index]=min(v(100:end));
Vpp=Vp-Vr;

%Derivative of AP
dVdt=zeros(length(v)-1,1);
for i=1:(length(v)-1)
    dVdt(i)=(v(i+1)-v(i))/(t(i+1)-t(i));
end

% Detect repolarization abnormalities
t_ind_repol= find(t>100,1); 
if sum(dVdt(t_ind_repol:end)>0.01)>1 || Vr > -70
    APDxx=NaN; Vpp=NaN; dVdt_max=NaN; Vr=NaN; Vp=NaN;
    EAD=1;
else
    [dVdt_max, dVdt_max_index]=max(dVdt);
    t_maxderv=t(dVdt_max_index);

    V_xxrep=Vp-(xx/100)*Vpp;
    v_rep=v(Vp_index:end); 
    [dif,vxx_index_cut]=min(abs(v_rep-V_xxrep));
    vxx_index=vxx_index_cut+(Vp_index-1); 
    t1=t(vxx_index); v1=v(vxx_index);
    if v(vxx_index)>V_xxrep
        t2=t(vxx_index+1); v2=v(vxx_index+1);
    end
    if v(vxx_index)<V_xxrep
        t2=t(vxx_index-1); v2=v(vxx_index-1);
    end
    t_vxx=t1+(V_xxrep-v1)*((t2-t1)/(v2-v1));
    APDxx=t_vxx-t_maxderv;
    EAD=0;
end
end

function [CTDxx,cai_syst,cai_diast, EAD]=CTD(t,cai,xx, flag)  
%Ca2+-transient duration at xx% of the initial base value (xx, eg.90%)-Passini et al.
%For 1 APD:
%%%% t:time vector
%%%% cai: vector of intracellular Ca2+
%%%% t0: stimulus time
%%%% flag indicates if an EAD has been detected previously

[cai_syst, cai_sys_index]=max(cai);
t_ind_rep=find(t>300,1); 
[cai_diast, cai_diast_index]=min(cai(t_ind_rep:end));

dCaidt=diff(cai)./diff(t);

% EADs
if sum(dCaidt(t_ind_rep:end)>0.0001)>1 | flag == 1
    CTDxx=NaN;
    cai_syst=NaN;
    cai_diast=NaN;
    EAD = 1;
else
    
    Acai=cai_syst-cai_diast;
    
    %CaTD
    cat_xxrec=cai_syst-(xx/100)*Acai;
    t_catxx = interp1(cai(cai_sys_index:end),t(cai_sys_index:end),cat_xxrec);
    [~,ind]=max(dCaidt);
    CTDxx=t_catxx-t(ind);
    EAD = 0;

end
end