function rr = simPAF_gen_AF_RR_intervals(lamba,nRR)
% rr = simPAF_gen_AF_RR_intervals() returns AF RR series modelled according
% to the paper by Corino et al. An atrioventricular node model for analysis
% of the ventricular response during atrial fibrillation. IEEE Transactions
% on Biomedical Engineering. 2011, 58(12), 3386-3395.
%
% Copyright (C) 2011 Valentina D. A. Corino*, Frida Sandberg**,
% Luca T. Mainardi*, Leif Sornmo**
% *Department of Bioengineering, Politecnico di Milano;
% **Department of Biomedical Engineering and
% Center of Integrative Electrocardiology, Lund University;
%
% Available under the GNU General Public License version 3
% - please see the accompanying file named "LICENSE"
%
if nRR < 70
time = 1;
else
time = nRR/70;
end
prob_alpha = 0.6;
P = [0.15 0.25];
difftau = 0.2;
slope_beta = 10;
tau1_in=polyval(P,1);
tau1=polyval(P,0:.001:3);
P2=[P(1) P(2)+difftau];
tau2=polyval(P2,0:.001:3);
tau2_in=polyval(P2,1);
Vt=-40;
Vr=-90;
dVdt=0;
aa=exprnd(1/lamba,1000000,1);
t_a=cumsum(aa);
aa=aa(t_a<time*60);
atr_t=0;
avj_tmA=0;
vtr_tmA=0;
phase='phase4';
Vm=Vr;
Ts=0.001;
t=0;
nextAA=aa(1);
i_R=0;
vtr_t=0;
vtr_nA=0;
avj_t=0;avj_nA=0;
R=zeros(length(aa),1);
nAA=0;
prob=zeros(length(aa),1);
prob1=round(length(aa)*prob_alpha);
prob(1:prob1)=1;
prob=prob(randperm(length(prob)));
prob_vera=zeros(size(prob));
tautau=zeros(size(aa));
tempo_x=0:.001:3;
beta=-slope_beta*tempo_x+1;
beta(beta<0)=0;
prob2=round(length(aa)*beta);
Nbeta=length(find(prob2>0));
prob_beta=zeros(length(aa),length(Nbeta));
for ii=1:Nbeta
prob_beta(1:prob2(ii),ii)=1;
prob_beta(:,ii)=prob_beta(randperm(length(aa)),ii);
end
blocked_beta=zeros(length(aa),1);
blocked_beta_t=zeros(length(tempo_x),1);
non_blocked_beta_t=zeros(length(tempo_x),1);
tempo_tot=zeros(length(aa),1);
while nAA<length(aa)
%%% function UpdateAtTs
% At every msec this is what happens: Vm is linearly increased
% by dVdt*Ts if in phase4; or the avj_t is increased
t=t+Ts; % update time
atr_t=atr_t+Ts;
vtr_t=vtr_t+Ts;
avj_tmA=avj_tmA+Ts;
vtr_tmA=vtr_tmA+Ts;
% update timers
if strcmp(phase,'phase4') % otherwise it's in refractory period
Vm=Vm+dVdt*Ts;
else
avj_t=avj_t+Ts;
end
%%% function AnteHitAVJ
% an AA arrives at the AVJ: Vm is increased of deltaV if in phase4 or
% tau (=RP) is prolonged
if atr_t>=nextAA
% t1=atr_t;
atr_t=0;
nAA=nAA+1; % increments AF counter
if nAA<length(aa)
nextAA=aa(nAA+1); % next AA
end
if strcmp(phase,'phase4') % otherwise it's in refractory period
beats=find(R>0);
if beats>0
tempo=t-R(beats(end))-tau;
[m, pos]=min(abs(tempo_x-tempo));
if pos<=Nbeta
if prob_beta(nAA,pos)==0
deltaV=(Vt-Vr)+1;
Vm=Vm+deltaV;
non_blocked_beta_t(pos)=non_blocked_beta_t(pos)+1;
else blocked_beta(nAA)=1;
blocked_beta_t(pos)=blocked_beta_t(pos)+1;
tempo_tot(nAA)=tempo;
end
else
deltaV=(Vt-Vr)+1;
Vm=Vm+deltaV;
non_blocked_beta_t(pos)=non_blocked_beta_t(pos)+1;
end
else
deltaV=(Vt-Vr)+1;
Vm=Vm+deltaV;
end
end
end
%%% function VtrSense
% there is a wave in the ventricle (vtr_nA>0) and the ventricle
% is not in refractory period (vtr_tmA>=AntDly)
if vtr_nA>0
i_R=i_R+1;
R(i_R)=t;
vtr_tmA=0;
vtr_t=0;
vtr_nA=0;
end
%%% function ActivateAVJ + StartAVJref if in phase4:
% when Vm>=Vt AVJ is activated (avj_nA=avj_nA+1;) and then the
% RP starts in AVJ
% OR StartAVJph4 if in phase0
% when avj_t>tau, i.e. the RP of AVJ finishes, phase4 again
if strcmp(phase,'phase4') % otherwise it's in refractory period
if Vm>=Vt
%%% ActivateAVJ
avj_nA=avj_nA+1;
vtr_nA=vtr_nA+1;
phase='phase0';
avj_t=0;
if prob(nAA)==0
beats=find(R>0);
if beats>0
new_RR=t-R(beats(end));
[m xx]=min(abs(tempo_x-new_RR));
tau=tau2(xx);
prob_vera(nAA)=2;
else
tau=tau2_in;
prob_vera(nAA)=2;
end
else
beats=find(R>0);
if beats>0
new_RR=t-R(beats(end));
[m, xx]=min(abs(tempo_x-new_RR));
tau=tau1(xx);
prob_vera(nAA)=1;
else
tau=tau1_in;
prob_vera(nAA)=1;
end
end
tautau(nAA)=tau;
end
elseif exist('tau1','var')==1
if avj_t>tau
%%% StartAVJph4
phase='phase4';
Vm=Vr;
end
end
end
R=R(R>0);
rr=diff(R);