function [biomark, EADsInds, numEADs] = calculateBiomarkers_calibration(s, s_Iks, s_Ikr, s_Ik1)
% Inputs:
%
% s: struc with the results of the simulation of the population in
% control conditions
% s_Iks: struc with the results of the simulation of the population with
% IKs*0.1
% s_Ikr: struc with the results of the simulation of the population with
% IKr*0.3
% s_Ik1: struc with the results of the simulation of the population with
% IKs*0.9
%
% biomark=[apd40 apd50 apd90 tri90_40 dvdtmax vpeak rmp...
%             ctd50 ctd90 cai_syst cai_diast nai_syst  delta_Iks delta_Ikr
%             delta_Ik1 EAD Na_alteration]; 
      

biomark=zeros(length (s.result),17,3);
for j=1:length (s.result)
    for i=1:size(s.result{1,1},1)
        t=s.result{j,1}{i,1};
        v=s.result{j,2}{i,1}(:,1);
        t_Iks=s_Iks.result{j,1}{i,1};
        v_Iks=s_Iks.result{j,2}{i,1}(:,1);
        t_Ikr=s_Ikr.result{j,1}{i,1};
        v_Ikr=s_Ikr.result{j,2}{i,1}(:,1);
        t_Ik1=s_Ik1.result{j,1}{i,1};
        v_Ik1=s_Ik1.result{j,2}{i,1}(:,1);
        cai=s.result{j,2}{i,1}(:,6)*10^3;% en uM
        nai=s.result{j,2}{i,1}(:,2);% en mM

        
        [apd90,~,dvdtmax,rmp,vpeak, EAD]=APD(t,v,90);
        [apd90_Iks,~,~,~,~, ~]=APD(t_Iks,v_Iks,90);
        [apd90_Ikr]=APD(t_Ikr,v_Ikr,90);
        [apd90_Ik1]=APD(t_Ik1,v_Ik1,90);
        delta_Iks=(apd90_Iks-apd90)/apd90*100;
        delta_Ikr=(apd90_Ikr-apd90)/apd90*100;
        delta_Ik1=(apd90_Ik1-apd90)/apd90*100;
        apd50=APD(t,v,50);
        apd40=APD(t,v,40);
        tri90_40=apd90-apd40;
        
        [ctd50,cai_syst,cai_diast]=CTD(t,cai,50, EAD);
        ctd90=CTD(t,cai,90, EAD);
        [nai_syst,nai_sys_index]=max(nai);
        if t(nai_sys_index) > 900
            Na_alteration=1;
        else
            Na_alteration=0;
        end
        %                1     2      3       4       5       6   7
        biomark(j,:,i)=[apd40 apd50 apd90 tri90_40 dvdtmax vpeak rmp...
            ctd50 ctd90 cai_syst cai_diast nai_syst  delta_Iks delta_Ikr delta_Ik1 EAD Na_alteration];
        %     8     9         10      11        12       13       14          15     16       17
    end
 [EADsinds1,~]=find(biomark(:,10,1)==1|biomark(:,11,1)==1);
 [EADsinds2,~]=find(biomark(:,10,2)==1|biomark(:,11,2)==1);
 [EADsinds3,~]=find(biomark(:,10,3)==1|biomark(:,11,3)==1);
 
 EADsInds=unique([EADsinds1; EADsinds2; EADsinds3]);
 numEADs=numel(EADsInds);
 
 end
end
%%
function [APDxx,Vpp,dVdt_max,Vr,Vp, EAD]=APD(t,v,xx)
%function [APDxx,t_vxx,t_maxderv]=APD(t,v,xx)    

%APD at different repolarization % (xx, eg 90%)
%Other parameters: amplitude (Vpp), max upstroke (dVdt_max), Vrest, Vpeak
    
[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
% EADs detection
t_ind_repol= find(t>150,1);
if sum(dVdt(t_ind_repol:end)>0.05)>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]=CTD(t,cai,xx, flag)  
%Ca2+-transient duration at xx% of the initial base value (xx, eg.90%)-Passini et al.
%%%% 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 detection
if sum(dCaidt(t_ind_rep:end)>0.0001)>1 || flag == 1
    CTDxx=NaN;
    cai_syst=NaN;
    cai_diast=NaN;
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);
end
end