/* afvp.c */
/* Copyright (c)2006 by Jie Lian, All Rights Reserved */
/* This program simulates the interaction between AF and Vp. */
/* See Lian J., Mussig D., Lang V. Computer Modeling of Ventricular */
/* Rhythm During Atrial Fibrillation and Ventricular Pacing. */
/* IEEE Transactions on Biomedical Engineering 53(8): 1512-1520, 2006. */
/* Contact J. Lian jie.lian AT biotronik DOT com) */
/* */
/* This program is free software; you can redistribute it and/or modify */
/* it under the terms of the GNU General Public License as published by */
/* the Free Software Foundation; either version 2 of the License, or */
/* (at your option) any later version. */
/* */
/* This program is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
/* GNU General Public License for more details. */
/* */
/* You should have received a copy of the GNU General Public License */
/* along with this program; if not, write to the Free Software Foundation */
/* Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
/* */
/* afvp.c and its dependents are freely available from Physionet - */
/* http://www.physionet.org/ - please report any bugs to the author above. */
/* Last modified: Jul 12th 2006 by Gari Clifford: */
/* - Adapted random number calls to make it ANSI C compliant */
/* - Added random number seeding and extra command line arguments */
/* - Removed BI (standby ventricular pacing basic interval) command line */
/* option */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h> // For random()
#include "afvp.h" // header file for the AF-VP model
#include "conio.h" // sometimes a missing system file
// definition of global model variables
EnvStruct env; // simulation environment
AtrStruct atr; // atrium model
AvjStruct avj; // AVJ model
VtrStruct vtr; // ventricle model
VeleStruct vele; // ventricular electrode model
unsigned int rseed; // seed for random number generator
// -can be passed as an argument
int main(int argc, char *argv[])
{
time_t utc;
if (argc<2){
fprintf(stdout,"Usage: %s <config_file> <random number seed>\n",argv[0]);
fprintf(stdout,"\n config file example: http://www.physionet.org/physiotools/afvp/config.txt \n");
// fprintf(stdout," BI = standby ventricular pacing basic interval in ms (default = 800 ms) \n");
fprintf(stdout," rseed = seed for random number generator - specify a long int for repeatability \n\n");
exit(1);
}
ReadConfig(argv[1]); // read model configuration file
//if (argc==3) vele.BI = atoi(argv[2])/1000.0; // overwrite BI
if (argc>2){
rseed = atol(argv[2]); // random seed
}
else {
/* use the system clock to give it a random number */
time(&utc);
rseed = (long int) utc;
}
if (argc>3) fprintf(stderr,"Ignoring input after 2nd argument\n");
ModelInit(); // initialize model parameters
// Main loop, keep running until MAX_RR beats or MAX_TIME timeout
while ((env.beat<env.MAX_RR) && (env.t<env.MAX_TIME)) {
UpdateAtTs(); // model update at sampling time
if (vtr.tmA[0]/vtr.AntDly+vtr.tmR[0]/vtr.RetDly>=1.0)
VtrFusion(); // detection of ventricular fusion
if (atr.t>=atr.NextAA)
AnteHitAvj(); // antegrade AF impluse hits AVJ
if (vtr.nR>0 && vtr.tmR[0]>=vtr.RetDly)
RetrHitAvj(); // VP-induced retrograde impulse hits AVJ
if (vtr.nA>0 && vtr.tmA[0]>=vtr.AntDly)
VtrSense(); // there is an ante conducted Vs
if (vele.t>=vele.NextVV)
VtrPace(); // timeout for a scheduled Vp
if (avj.tmA[0][0]>=avj.tmA[0][1])
AnteEscAvj(); // ante AV delay timeout and invade vtr
if (avj.tmR[0][0]>=avj.tmR[0][1])
RetrEscAvj(); // retr AV delay timeout and invade atr
if (avj.phase==PHASE4) { // check if AVJ depol & starts ref
if (avj.Vm>=avj.Vt) {
ActivateAvj();
StartAvjRef();
}
}
else { // check if AVJ ref ends & return to phase 4
if (avj.t>=avj.ref)
StartAvjPh4();
}
}
ModelExit();
} /*** end main ***/
/*********************************************************************
ErrMsg: display error message and abort simulation
Input: msg - message to be displayed
errcode - exit error code
*********************************************************************/
void ErrMsg(char *msg, int errcode)
{
printf("%s\n",msg);
exit(errcode);
}
/*********************************************************************
LogEvent: log event time and type
Input: fp - log file pointer
msg - message to be logged
*********************************************************************/
void LogEvent(FILE *fp, char *msg)
{
fprintf(fp,"T=%.3lf '%s ",env.t,msg);
if (avj.phase==PHASE4) fprintf(fp,"@PH4' ");
else fprintf(fp,"@REF' ");
}
/*********************************************************************
LogStat: log statistics after simulation
Input: fp - log file pointer
*********************************************************************/
void LogStat(FILE *fp)
{
long sum;
fprintf(fp,"\n----------------------------------------------------\n");
fprintf(fp,"Random seed = %d\n",rseed);
fprintf(fp,"Total time (s): %.3lf\n", env.t);
fprintf(fp,"Total beat: %ld\n", env.beat);
fprintf(fp,"Vs Count = %ld [%.2lf%%]\n",
vele.cnt[SENSE],100.0*vele.cnt[SENSE]/env.beat);
fprintf(fp,"Vp Count = %ld [%.2lf%%]\n",
vele.cnt[PACE], 100.0*vele.cnt[PACE]/env.beat);
sum = atr.nAf[REFRACTORY]+atr.nAf[PHASE4];
fprintf(fp,"Af Count = %ld(ref)+%ld(ph4) = %ld\n",
atr.nAf[REFRACTORY], atr.nAf[PHASE4], sum);
fprintf(fp,"Mean RR interval = %.3lf <%dppm>\n",
env.t/env.beat, (int)(env.beat*60/env.t));
fprintf(fp,"Mean PP interval = %.3lf <%dppm>\n",
env.t/sum, (int)(sum*60/env.t));
sum = vtr.nR+avj.nR;
fprintf(fp,"Active retr waves = %ld [%.2lf%%]\n",
sum, 100.0*sum/env.beat);
sum = vtr.nF[REFRACTORY]+vtr.nF[PHASE4];
fprintf(fp,"Vtr fusions = %ld(ref)+%ld(ph4) = %ld [%.2lf%%]\n",
vtr.nF[REFRACTORY], vtr.nF[PHASE4], sum, 100.0*sum/env.beat);
sum = avj.nF[REFRACTORY]+avj.nF[PHASE4];
fprintf(fp,"AVJ fusions = %ld(ref)+%ld(ph4) = %ld [%.2lf%%]\n",
avj.nF[REFRACTORY], avj.nF[PHASE4], sum, 100.0*sum/env.beat);
sum = atr.nF[REFRACTORY]+atr.nF[PHASE4];
fprintf(fp,"Atr fusions = %ld(ref)+%ld(ph4) = %ld [%.2lf%%]\n",
atr.nF[REFRACTORY], atr.nF[PHASE4], sum, 100.0*sum/env.beat);
sum = avj.nB[RETROGRADE]+vtr.nB[RETROGRADE];
fprintf(fp,"Retrgrade Block = %ld(avj)+%ld(vtr) = %ld [%.2lf%%]\n",
avj.nB[RETROGRADE], vtr.nB[RETROGRADE], sum, 100.0*sum/env.beat);
}
/*********************************************************************
ReadConfig: read configuration file
Input: fn - config filename
*********************************************************************/
void ReadConfig(char *fn)
{
double d;
FILE *fp;
char str[BUFLIM];
if ((fp=fopen(fn,"rt"))==NULL) ErrMsg("Cannot open config file",1);
// read parameters of the simulation environment
ReadLine(fp,-1,env.fnRR);
ReadLine(fp,-1,env.fnAA);
ReadLine(fp,-1,env.fnAV);
ReadLine(fp,-1,env.fnLOG);
env.MAX_RR = (long)ReadLine(fp,500,str);
env.MAX_TIME = ReadLine(fp,1000,str);
env.Ts = ReadLine(fp,0.001,str);
env.RR0 = ReadLine(fp,1.0,str);
// read parameters of the atrium model
atr.AA_MODEL = (int)ReadLine(fp,0,str);
atr.lambda = ReadLine(fp,10,str);
atr.AAstd = ReadLine(fp,0,str);
atr.dVmean = ReadLine(fp,33,str);
atr.dVstd = ReadLine(fp,0,str);
atr.AtrDly = ReadLine(fp,0.03,str);
atr.S1S2 = ReadLine(fp,0.4,str);
atr.S2S3 = ReadLine(fp,0.4,str);
// read parameters of the AVJ model
avj.Vt = ReadLine(fp,-40,str);
avj.Vr = ReadLine(fp,-90,str);
avj.dVdt = ReadLine(fp,30,str);
avj.MinAVDa = ReadLine(fp,0.07,str);
avj.MinAVDr = ReadLine(fp,0.07,str);
avj.alpha = ReadLine(fp,0.13,str);
avj.tau_c = ReadLine(fp,0.1,str);
avj.MinRef = ReadLine(fp,0.05,str);
avj.beta = ReadLine(fp,0.25,str);
avj.tau_r = ReadLine(fp,0.5,str);
avj.Ref_std = ReadLine(fp,0.0,str);
avj.delta = ReadLine(fp,10,str);
avj.theta = ReadLine(fp,10,str);
// read parameters of the ventricle model
vtr.AntDly = ReadLine(fp,0.05,str);
vtr.RetDly = ReadLine(fp,0.15,str);
vtr.ref = ReadLine(fp,0.1,str);
// read parameters of the RV electrode model
vele.VP_MODEL = (int)ReadLine(fp,0,str);
vele.BI = ReadLine(fp,1.0,str);
fclose(fp);
LogConfig(stdout); // display model parameters on terminal
}
/*********************************************************************
ReadLine: read a line from the config file
Input: fp - pointer to the config file
v - default value
str - char pointer (string read after '=' sign)
Output: 'val' read from the line if found, otherwise return 'v'
Note: Comment line starts with '%' sign
Value is separated from text by token '='
*********************************************************************/
double ReadLine(FILE *fp, double v, char *str)
{
char buf[BUFLIM], *p;
char tok[] = "="; // token followed by 'val'
double val; // value read from the line
*str = '\0'; // default: null string
do {
p = fgets(buf,sizeof(buf),fp); // read a line
if (p==NULL) return v; // force return on EOF
} while (strlen(buf)<=1 || buf[0]=='%'); // skip comment or blank lines
p = strstr(buf,tok); // search for "=" sign
if (p!=NULL) { // if found "="
p += strlen(tok); // skip "=" sign
sscanf(p,"%s",str); // read string after "=" sign
if (sscanf(p,"%lf",&val)>0) return val;// read value after "=" sign
} // otherwise
return v; // return default value
}
/*********************************************************************
LogConfig: log model configuration
Input: fp - output file pointer
*********************************************************************/
void LogConfig(FILE *fp)
{
fprintf(fp,"----------------------------------------------------\n");
fprintf(fp,"*** Simulation Environment ***\n");
fprintf(fp,"fnRR=%s\nfnAA=%s\nfnLog=%s\nfnAV=%s\n",
env.fnRR,env.fnAA,env.fnAV,env.fnLOG);
fprintf(fp,"MAX_RR=%ld\nMAX_TIME=%lf\nTs=%lf\nRR0=%lf\n",
env.MAX_RR,env.MAX_TIME,env.Ts,env.RR0);
fprintf(fp,"\n*** Atrium Model ***\n");
fprintf(fp,"AA_MODEL=%d\nlambda=%lf\nAAstd=%lf\n",
atr.AA_MODEL,atr.lambda,atr.AAstd);
fprintf(fp,"dVmean=%lf\ndVstd=%lf\nAtrDly=%lf\nS1S2=%lf\nS2S3=%lf\n",
atr.dVmean,atr.dVstd,atr.AtrDly,atr.S1S2,atr.S2S3);
fprintf(fp,"\n*** AVJ Model ***\n");
fprintf(fp,"Vt=%lf\nVr=%lf\ndVdt=%lf\nMinAVDa=%lf\nMinAVDr=%lf\n",
avj.Vt,avj.Vr,avj.dVdt,avj.MinAVDa,avj.MinAVDr);
fprintf(fp,"alpha=%lf\ntau_c=%lf\nMinRef=%lf\nbeta=%lf\n",
avj.alpha,avj.tau_c,avj.MinRef,avj.beta);
fprintf(fp,"tau_r=%lf\nRef_std=%lf\ndelta=%lf\ntheta=%lf\n",
avj.tau_r,avj.Ref_std,avj.delta,avj.theta);
fprintf(fp,"\n*** Ventricle Model ***\n");
fprintf(fp,"AntDly=%lf\nRetDly=%lf\nref=%lf\n",
vtr.AntDly,vtr.RetDly,vtr.ref);
fprintf(fp,"\n*** Ventricular Electrode Model ***\n");
fprintf(fp,"VP_MODEL=%d\nBI=%lf\n",vele.VP_MODEL,vele.BI);
fprintf(fp,"----------------------------------------------------\n\n");
}
/*********************************************************************
ModelInit: initialization of model parameters
*********************************************************************/
void ModelInit()
{
int i;
if ((env.fprr=fopen(env.fnRR,"wt"))==NULL) // output file of RR intervals
ErrMsg("Cannot open output file of RR intervals",1);
if ((env.fpaa=fopen(env.fnAA,"wt"))==NULL) // output file of AA intervals
ErrMsg("Cannot open output file of AA intervals",1);
if ((env.fpav=fopen(env.fnAV,"wt"))==NULL) // output file of AVJ status
ErrMsg("Cannot open output file of AVJ status",1);
if ((env.fplog=fopen(env.fnLOG,"wt"))==NULL)// output file of event log
ErrMsg("Cannot open output file of event log",1);
env.beat = 0; // clear beat counter
env.t = 0; // reset simulation time clock
avj.phase = PHASE4; // set AVJ phase = PHASE4
avj.Vm = avj.Vr; // set Vm to Vr
avj.V4 = avj.dVdt * env.Ts; // initialize V4 (mV/ms)
avj.RT0 = 0; // AVJ start recovery
avj.nA = avj.nR = 0; // no activation wave in AVJ
vtr.nA = vtr.nR = 0; // no activation wave in vtr
for (i=0;i<MAX_AVJ_TM;i++) { // disable AVJ activation timers
avj.tmA[i][0] = -1;
avj.tmR[i][0] = -1;
}
for (i=0;i<MAX_VENT_TM;i++) { // disable vtr activation timers
vtr.tmA[i] = -1;
vtr.tmR[i] = -1;
}
vele.cnt[PACE] = 0; // clear Vp counter
vele.cnt[SENSE] = 0; // clear Vs counter
atr.nAf[0] = atr.nAf[1] = 0; // clear AF event counters
atr.nF[0] = atr.nF[1] = 0; // clear atrial fusion counters
avj.nF[0] = avj.nF[1] = 0; // clear AVJ fusion counters
vtr.nF[0] = vtr.nF[1] = 0; // clear ventricle fusion counters
vtr.tmA0 = 0; // clear last ante AV timeout time
vtr.tmR0 = 0; // clear last Vp delivery time
// Get the first ventricular pacing interval
vele.NextVV = GetVpIntv(env.RR0,vele.BI,vele.BI,vele.VP_MODEL);
// randomize(); // Initialize randomization
srandom(rseed); // Initialize randomization with seeded random number
// Get the first random interval for next AF depol.
atr.NextAA = GetAfIntv(atr.lambda,atr.AAstd,atr.AA_MODEL);
fprintf(env.fpaa,"%lf\n",atr.NextAA);
LogConfig(env.fplog); // log config parameters
}
/*********************************************************************
ModelExit: exit simulation
*********************************************************************/
void ModelExit()
{
LogStat(stdout); // output statistics to stdout
LogStat(env.fplog); // output statistics to log file
fclose(env.fprr);
fclose(env.fpaa);
fclose(env.fpav);
fclose(env.fplog);
}
/*********************************************************************
UpdateAtTs: model update at sampling time
*********************************************************************/
void UpdateAtTs()
{
int i;
env.t += env.Ts; // update simulation timer
atr.t += env.Ts; // update atrial timer
vtr.t += env.Ts; // update ventricle timer
vele.t += env.Ts; // update Vp timer
// update AVJ and ventricle activation timers
for (i=0;i<avj.nA;i++) avj.tmA[i][0] += env.Ts;
for (i=0;i<avj.nR;i++) avj.tmR[i][0] += env.Ts;
for (i=0;i<vtr.nA;i++) vtr.tmA[i] += env.Ts;
for (i=0;i<vtr.nR;i++) vtr.tmR[i] += env.Ts;
if (avj.phase==PHASE4) { // during phase 4
avj.Vm += avj.V4; // Vm linear rise from Vr to Vt
if (avj.Vm>=avj.Vt) avj.dir = ANTEGRADE;
}
else // during refractory period
avj.t += env.Ts; // update AVJ refractory timer
}
/*********************************************************************
AnteHitAvj: event service when antegrade impulse hits AVJ
*********************************************************************/
void AnteHitAvj()
{
double ext;
long r;
LogEvent(env.fplog,"AnteHitAvj");
atr.t = 0; // reset atrium timer
atr.nAf[avj.phase]++; // increment AF impulse counter
// Get the random interval for next AF depolarization
atr.NextAA = GetAfIntv(atr.lambda,atr.AAstd,atr.AA_MODEL);
fprintf(env.fpaa,"%lf\n",atr.NextAA);
if (avj.phase==PHASE4) { // AVJ in phase IV
//srandom(-1000);
r = (long) random(); // simulate variation of dV
atr.dV = gaussian(atr.dVmean,atr.dVstd,ran2(&r),ran2(&r));
if (atr.dV>0) avj.Vm += atr.dV; // delta increment of Vm
if (avj.Vm>=avj.Vt) avj.dir = ANTEGRADE;
fprintf(env.fplog, // log after AnteHitAvj() service @PH4
"NextAA=%.3lf, Vm=%.2lf+%.2lf=%.2lf\n",
atr.NextAA, avj.Vm-atr.dV, atr.dV, avj.Vm);
}
else { // AVJ in refractory period
r = atr.dV/(avj.Vt-avj.Vr); // relative strength of the AF impulse
if (r>1) r = 1.0; // with upper limit 1.0
r = pow(avj.t/avj.ref,avj.theta)*pow(r,avj.delta);
ext = avj.MinRef*r;
avj.ref += ext; // ref extension due to concealed AF
fprintf(env.fplog, // log after AnteHitAvj() service @REF
"NextAA=%.3lf, ref=%.3lf+%.3lf=%.3lf\n",
atr.NextAA, avj.ref-ext, ext, avj.ref);
}
}
/*********************************************************************
RetrHitAvj: event service when retrograde impulse hits AVJ
*********************************************************************/
void RetrHitAvj()
{
int i;
double r,ext;
LogEvent(env.fplog,"RetrHitAvj");
StopVtrTm(RETROGRADE); // stop a retr vtr activation timer
if (avj.phase==PHASE4) { // retr invasion in phase 4
if (avj.Vm<avj.Vt) { // retr AVJ depolarization
avj.Vm = avj.Vt; // suprathreshold AVJ immediately
avj.dir = RETROGRADE; // retrograde activation
fprintf(env.fplog, // log after RetHitAvj() service @PH4
"vtr.nR=%d, Vm=%.2lf Retr AVJ activation\n", vtr.nR, avj.Vm);
}
else { // simultaneous ante/retr depol (rare), no wave in AVJ
avj.dir = BIDIRECTION;
fprintf(env.fplog, // log after RetHitAvj() service @PH4
"vtr.nR=%d, Vm=%.2lf Bidir AVJ activation\n", vtr.nR, avj.Vm);
}
}
else { // retr invasion in refractory period
if (avj.nA>0) { // presence of ante wave in AVJ
StopAvjTm(ANTEGRADE); // stop the ante AVJ timer
avj.nF[avj.phase]++; // incr AVJ fusion counter
}
avj.nB[RETROGRADE]++; // incr counter of blocked retr waves
r = pow(avj.t/avj.ref,avj.theta); // strong impulse
ext = avj.MinRef*r;
avj.ref += ext; // ref extension by concealed impulse
fprintf(env.fplog, // log after RetHitAvj() service @REF
"vtr.nR=%d, avj.nA=%d, nB=%d, nF=%d, ref=%.2lf+%.2lf=%.2lf\n",
vtr.nR,avj.nA,avj.nB[RETROGRADE],avj.nF[REFRACTORY],
avj.ref-ext,ext,avj.ref);
}
}
/*********************************************************************
ActivateAvj: AVJ depolarization service
*********************************************************************/
void ActivateAvj()
{
LogEvent(env.fplog,"ActivateAvj");
avj.RT = env.t-avj.RT0; // calculated AVJ recovery time
avj.AVDa = avj.MinAVDa+avj.alpha*exp(-avj.RT/avj.tau_c); // ante AV delay
avj.AVDr = avj.MinAVDr+avj.alpha*exp(-avj.RT/avj.tau_c); // retr AV delay
if (avj.dir==ANTEGRADE) { // antegrade AVJ activation
if (avj.nR==0) { // no existing retr activation in AVJ
StartAvjTm(ANTEGRADE); // start a new timer for ante AV delay
fprintf(env.fpav, // record AVJ status
"%.3lf %d %.3lf %.3lf ", env.t, avj.dir, avj.RT, avj.AVDa);
fprintf(env.fplog, // log after ActivateAvj() service
"start AVDa = %.3lf, avj.nA=%d\n", avj.AVDa, avj.nA);
}
else { // collide w/ existing retr wave
StopAvjTm(RETROGRADE); // stop a retr AVJ activation timer
avj.nF[avj.phase]++; // incr AVJ collision counter (Phase 4)
fprintf(env.fpav, // record AVJ status
"%.3lf %d %.3lf %d ", env.t, avj.dir, avj.RT, -1);
fprintf(env.fplog, // log after ActivateAvj() service
"collide w/ retr wave, avj.nF(%d)=%ld\n", avj.phase, avj.nF[avj.phase]);
}
}
else if (avj.dir==RETROGRADE) { // retrograde AVJ activation
if (avj.nA==0) { // no existing ante activation in AVJ
StartAvjTm(RETROGRADE); // start a new timer for retr AV delay
fprintf(env.fpav, // record AVJ status
"%.3lf %d %.3lf %.3lf ", env.t, avj.dir, avj.RT, avj.AVDr);
fprintf(env.fplog, // log after ActivateAvj() service
"start AVDr = %.3lf, avj.nR=%d\n", avj.AVDr, avj.nR);
}
else { // collide w/ existing ante wave in AVJ
StopAvjTm(ANTEGRADE); // stop an ante AVJ activation timer
avj.nF[avj.phase]++; // incr AVJ collision counter (phase 4)
fprintf(env.fpav, // record AVJ status
"%.3lf %d %.3lf %d ", env.t, avj.dir, avj.RT, -1);
fprintf(env.fplog, // log after ActivateAvj() service
"collide w/ ante wave, avj.nF(%d)=%ld\n", avj.phase, avj.nF[avj.phase]);
}
}
else { // simultaneous bidirectional depol
avj.nF[avj.phase]++; // incr AVJ collision counter (phase 4)
fprintf(env.fpav, // record AVJ status
"%.3lf %d %.3lf %d ", env.t, avj.dir, avj.RT, -1);
fprintf(env.fplog, // log after ActivateAvj() service
"bidir depolarization, avj.nF(%d)=%ld\n", avj.phase, avj.nF[avj.phase]);
}
}
/*********************************************************************
StartAvjRef: AVJ starts refractory period
*********************************************************************/
void StartAvjRef()
{
long r;
// calculate AVJ refractory period
avj.ref = avj.MinRef+avj.beta*(1-exp(-avj.RT/avj.tau_r));
// Assume AVJ refractory period does not exceed the AVD
if (avj.nA>0 && avj.ref<avj.AVDa) avj.ref = avj.AVDa; // limited by AVD
if (avj.nR>0 && avj.ref<avj.AVDr) avj.ref = avj.AVDr; // limited by AVD
//srandom(-1000);
r = (long) random();
avj.ref = gaussian(avj.ref,avj.Ref_std,ran2(&r),ran2(&r));
avj.phase = REFRACTORY; // update AVJ phase to REFRACTORY
avj.t = 0; // start AVJ refractory timer
fprintf(env.fpav,"%.3lf\n", avj.ref); // record AVJ status
fprintf(env.fplog, // log after StartAvjRef() service
"T=%.3lf 'Start AVJ REF' RT=%.3lf, ref=%.3lf\n",
env.t, avj.RT, avj.ref);
}
/*********************************************************************
StartAvjPh4: AVJ starts phase 4
*********************************************************************/
void StartAvjPh4()
{
avj.phase = PHASE4; // update AVJ phase to PHASE4
avj.Vm = avj.Vr; // set Vm to Vr
avj.RT0 = env.t; // AVJ starts recovery
fprintf(env.fplog, // log after StartAvjPh4() service
"T=%.3lf 'Start AVJ PH4' Vm=%.2lf, avj.nA=%d, avj.nR=%d\n",
env.t, avj.Vm, avj.nA, avj.nR);
}
/*********************************************************************
StartAvjTm: start a new activation timer in AVJ
Input: dir - direction of the activation
*********************************************************************/
void StartAvjTm(int dir)
{
if (dir==ANTEGRADE) { // a new timer for ante AV delay
avj.tmA[avj.nA][0] = 0; // start the timer
avj.tmA[avj.nA][1] = avj.AVDa; // set ante AV delay time
avj.nA++; // allow multiple ante waves in AVJ
}
else { // a new timer for retr AV delay
avj.tmR[avj.nR][0] = 0; // start the timer
avj.tmR[avj.nR][1] = avj.AVDr; // set retr AV delay time
avj.nR++; // allow multiple retr waves in AVJ
}
}
/*********************************************************************
StopAvjTm: stop an activation timer in AVJ
Input: dir - direction of the activation
*********************************************************************/
void StopAvjTm(int dir)
{
int i;
if (dir==ANTEGRADE) { // stop an ante AV delay timer
for (i=0;i<avj.nA-1;i++) { // update AVJ ante activation queue
avj.tmA[i][0] = avj.tmA[i+1][0]; // activation timer
avj.tmA[i][1] = avj.tmA[i+1][1]; // conduction delay
}
avj.nA--; // decr ante AVJ activation counter
avj.tmA[avj.nA][0] = -1; // disable an ante AVJ timer
}
else { // stop a retr AV delay timer
for (i=0;i<avj.nR-1;i++) { // update AVJ retr activation queue
avj.tmR[i][0] = avj.tmR[i+1][0]; // activation timer
avj.tmR[i][1] = avj.tmR[i+1][1]; // conduction delay
}
avj.nR--; // decr retro AVJ activation counter
avj.tmR[avj.nR][0] = -1; // disable a retr AVJ timer
}
}
/*********************************************************************
StartVtrTm: start a new activation timer in ventricle
Input: dir - direction of the activation
*********************************************************************/
void StartVtrTm(int dir)
{
if (dir==ANTEGRADE) { // a new timer for ante vtr delay
vtr.tmA[vtr.nA] = 0; // starts an ante activation in vtr
vtr.nA++; // incr ante vtr activation counter
if (vtr.nA>MAX_VENT_TM) ErrMsg("Too many ante waves in ventricle!",2);
}
else { // a new timer for retr vtr delay
vtr.tmR[vtr.nR] = 0; // starts a retr activation in vtr
vtr.nR++; // incr retr vtr activation counter
if (vtr.nR>MAX_VENT_TM) ErrMsg("Too many retr waves in ventricle!",2);
}
}
/*********************************************************************
StopVtrTm: stop an activation timer in ventricle
Input: dir - direction of the activation
*********************************************************************/
void StopVtrTm(int dir)
{
int i;
if (dir==ANTEGRADE) { // stop an ante vtr delay timer
for (i=0;i<vtr.nA-1;i++) // update ante vtr activation queue
vtr.tmA[i] = vtr.tmA[i+1];
vtr.nA--; // decr number of ante waves in vtr
vtr.tmA[vtr.nA] = -1; // stop an ante vtr activation timer
}
else { // stop a retr vtr delay timer
for (i=0;i<vtr.nR-1;i++) // update vtr retr activation queue
vtr.tmR[i] = vtr.tmR[i+1];
vtr.nR--; // decr retr vtr activation counter
vtr.tmR[vtr.nR] = -1; // stop a retr vtr activation timer
}
}
/*********************************************************************
AnteEscAvj: event service when ante wave escapes AVJ and invade vtr
*********************************************************************/
void AnteEscAvj()
{
LogEvent(env.fplog,"AnteEscAvj");
StopAvjTm(ANTEGRADE); // stop an ante AVJ activation timer
if (avj.nA==0) avj.RT0 = env.t; // AVJ starts recovery
if (env.t-vtr.tmA0<=vtr.ref) { // non-capture because vtr refr
vtr.nB[ANTEGRADE]++; // will not generate ante wave in vtr
fprintf(env.fplog, // log after AnteEscAvj() service
"avj.nA=%d, vtr.nB(a)=%d(V-Blk)\n", avj.nA, vtr.nB[ANTEGRADE]);
return; // force return
} // otherwise, capture vtr and ...
StartVtrTm(ANTEGRADE); // starts an ante activation in vtr
vtr.tmA0 = env.t; // update last ante captured time
fprintf(env.fplog, // log after AnteEscAvj() service
"avj.nA=%d, vtr.nA=%d\n", avj.nA, vtr.nA);
}
/*********************************************************************
RetrEscAvj: event service when retr wave escapes AVJ and invade atrium
*********************************************************************/
void RetrEscAvj()
{
double t;
LogEvent(env.fplog,"RetrEscAvj");
atr.nF[avj.phase]++; // incr atrial fusion counter
if (atr.NextAA-atr.t > atr.AtrDly) { // invade SA node before next impulse
t = GetAfIntv(atr.lambda,atr.AAstd,atr.AA_MODEL); // reset SA node
fprintf(env.fpaa,"%lf\n",t);
atr.NextAA += t; // update arrival time of next AF
} // otherwise, simply atrial fusion
StopAvjTm(RETROGRADE); // stop a retr AVJ activation timer
if (avj.nR==0) avj.RT0 = env.t; // AVJ starts recovery
fprintf(env.fplog, // log after RetrEscAvj() service
"avj.nR=%d, atr.nF(%d)=%ld, nextAA(%.3lf)+%.3lf=%.3lf\n",
avj.nR, avj.phase, atr.nF[avj.phase], atr.NextAA-t, t, atr.NextAA);
}
/*********************************************************************
VtrFusion: event service when detecting ventricular fusion
*********************************************************************/
void VtrFusion()
{
int i;
LogEvent(env.fplog,"VtrFusion");
vtr.nF[avj.phase]++; // incr ventricular fusion counter
StopVtrTm(ANTEGRADE); // stop an ante vtr activation timer
StopVtrTm(RETROGRADE); // stop a retr vtr activation timer
fprintf(env.fplog, // log after VtrFusion() service
"vtr.nF(%d)=%ld, vtr.nA=%d, vtr.nR=%d\n",
avj.phase, vtr.nF[avj.phase], vtr.nA, vtr.nR);
}
/*********************************************************************
RegEventV: register ventricular event (Vs/Vp) and update parameters
Input: event - PACE or SENSE
*********************************************************************/
void RegEventV(int event)
{
env.beat++; // increment beat counter
vele.VpFlg = event; // flag for sensed/paced beat
vele.cnt[event]++; // increment Vs/Vp counter
vele.RR = vtr.t; // record RR interval ended with Vs/Vp
// Get the next ventricular pacing interval
vele.NextVV = GetVpIntv(vele.RR,vele.BI,vele.NextVV,vele.VP_MODEL);
vele.t = 0; // reset Vp clock timer (inhibit Vp)
vtr.t = 0; // clear vtr time to measure RR interval
fprintf(env.fprr, // output RR interval and Vp/Vs flag
"%.3lf %d\n",vele.RR,event);
}
/*********************************************************************
VtrSense: event service for ventricular sense
*********************************************************************/
void VtrSense()
{
int i;
LogEvent(env.fplog,"VtrSense");
RegEventV(SENSE); // update parameters on Vs event
for (i=0;i<vtr.nA-1;i++) // update ante vtr activation timers
vtr.tmA[i] = vtr.tmA[i+1];
vtr.nA--;
vtr.tmA[vtr.nA] = -1; // stop an ante vtr activation timer
fprintf(env.fplog, // log after VtrSense() service
"vtr.nA=%d, nextVV=%.3lf, Vs#%ld, Beat#%ld, RR=%.3lf\n",
vtr.nA, vele.NextVV, vele.cnt[SENSE], env.beat, vele.RR);
}
/*********************************************************************
VtrPace: event service for ventricular pace
*********************************************************************/
void VtrPace()
{
LogEvent(env.fplog,"VtrPace");
RegEventV(PACE); // update parameters on Vp event
if (env.t-vtr.tmR0<=vtr.ref) { // non-capture because vtr refr
vtr.nB[RETROGRADE]++; // will not generate a retr wave in vtr
return; // force return
} // otherwise, Vp capture vtr and ...
StartVtrTm(RETROGRADE); // starts a retr activation in vtr
vtr.tmR0 = env.t; // update last captured Vp time
fprintf(env.fplog, // log after VtrPace() service
" vtr.nR=%d, nextVV=%.3lf, Vp#%ld, Beat#%ld, RR=%.3lf\n",
vtr.nR, vele.NextVV, vele.cnt[PACE], env.beat, vele.RR);
}
/*********************************************************************
GetAfIntv: get next Afib interval
*********************************************************************/
double GetAfIntv(double lambda, double sd, int method)
{
long init;
long k;
double next_afib_int, mean, next_afib_rate;
k = atr.nAf[0]+atr.nAf[1];
mean = 1/lambda;
//srandom(-1000);
init = (long) random();
switch (method) {
case 0: // random value from truncated exponential distribution
// next_afib_int = expdev(mean,&init);
next_afib_int = 0.05+expdev(mean-0.05,&init); // limit AA >= 50ms
break;
case 1: // inverse of random value from Poisson distribution
next_afib_rate = poidev(lambda,&init);
next_afib_rate++; // Ensure next_afib_rate>0 (avoid divided by 0)
next_afib_int = 1.0/next_afib_rate;
break;
case 2: // random value from uniform distribution
next_afib_int = mean+sd*mean*ran2(&init);
break;
case 3: // random value from Gaussian distribution
next_afib_int = gaussian(mean,sd*mean,ran2(&init),ran2(&init));
break;
case 4: // atrial pacing protocl with programmed S1S2
if (k < env.MAX_RR/2)
next_afib_int = 1.0/lambda;
else
next_afib_int = atr.S1S2;
break;
case 5: // atrial pacing protocl with programmed S1S2/S2S3
if (k < env.MAX_RR/2)
next_afib_int = 1.0/lambda;
else {
if (k%2==0) // even beat
next_afib_int = atr.S1S2;
else // odd beat
next_afib_int = atr.S2S3;
}
break;
default: // fixed rate, atrial flutter protocol
next_afib_int = 1.0/lambda;
break;
}
if (next_afib_int < 0.05) next_afib_int = 0.05; // ensure interval >= 50ms
return next_afib_int;
}