openCARP
Doxygen code documentation for the open cardiac electrophysiology simulator openCARP
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ap_analyzer.cc
Go to the documentation of this file.
1 // ----------------------------------------------------------------------------
2 // openCARP is an open cardiac electrophysiology simulator.
3 //
4 // Copyright (C) 2020 openCARP project
5 //
6 // This program is licensed under the openCARP Academic Public License (APL)
7 // v1.0: You can use and redistribute it and/or modify it in non-commercial
8 // academic environments under the terms of APL as published by the openCARP
9 // project v1.0, or (at your option) any later version. Commercial use requires
10 // a commercial license (info@opencarp.org).
11 //
12 // This program is distributed without any warranty; see the openCARP APL for
13 // more details.
14 //
15 // You should have received a copy of the openCARP APL along with this program
16 // and can find it online: http://www.opencarp.org/license
17 // ----------------------------------------------------------------------------
18 
34 #include <stdlib.h>
35 #include <string.h>
36 #include <math.h>
37 
38 #include "ap_analyzer.h"
39 #include "basics.h"
40 
41 namespace limpet {
42 
46 using ::opencarp::timer_manager;
47 
48 void shift_vm_trace(double vm, action_potential *AP, timer_manager *tm);
49 bool check_threshold(ap_event *e, timer_manager *tm);
50 bool check_mx_rate(ap_event *e, timer_manager *tm);
51 bool check_mx_val (ap_event *e, timer_manager *tm);
54 
55 
59 
60 
70 int
72 {
73  memset(AP,0,sizeof(action_potential));
74  AP->beat = -1;
75  // define which beat to use for calibration
77  AP->pmat = false;
78  int num_act_events = nActEvents;
79  int num_repol_events = nRepolEvents;
80  int num_ampl_events = nAmplEvents;
81 
82 #define BOGUS 9999
83  // analyze action potentials
84  // find activation marker
85  AP->acts.num_events = num_act_events,
86  AP->acts.events = (ap_event *)calloc(num_act_events,sizeof(ap_event));
87 
88  // set event type
89  ap_event *e = AP->acts.events;
90  e->c.m = X_THRESHOLD;
91  e->c.x_thresh.name = dupstr("LAT_xAPD90");
93  e->c.x_thresh.p_slope = true;
94  // initialize event values
95  e->s = BOGUS;
96  e->t = -10.0;
97  e->trace = AP->vm_trc; // either vm or dvm trace
98 
99  e++;
100  e->c.m = MX_VAL;
101  e->c.x_thresh.name = dupstr("LAT_dVdt");
102  e->c.x_thresh.p_slope = true;
103  // initialize event values
104  e->s = -1e10;
105  e->t = -10.;
106  e->trace = AP->dvm_trc; // find max of dvm trace
107 
108 
109  // find repolarisation marker
110  AP->repols.num_events = num_repol_events,
111  AP->repols.events = (ap_event *)calloc(num_repol_events,sizeof(ap_event));
112  // set event type
113  e = AP->repols.events;
114  e->c.m = X_THRESHOLD;
115  e->c.x_thresh.name = dupstr("REPOL_90");
117  e->c.x_thresh.p_slope = false;
118  // initialize event values
119  e->s = BOGUS;
120  e->t = -10.0;
121  e->trace = AP->vm_trc;
122 
123  e++;
124  e->c.m = X_THRESHOLD;
125  e->c.x_thresh.name = dupstr("REPOL_30");
127  e->c.x_thresh.p_slope = false;
128  // initialize event values
129  e->s = BOGUS;
130  e->t = -10.0;
131  e->trace = AP->vm_trc;
132 
133  // find extremal signal values
134  AP->ampls.num_events = num_ampl_events,
135  AP->ampls.events = (ap_event *)calloc(num_ampl_events,sizeof(ap_event));
136  // set event type
137  e = AP->ampls.events;
138  e->c.m = MX_VAL;
139  e->c.mx_val.name = dupstr("VM_MX");
140  e->c.mx_val.p_slope = true;
141  // initialize event values
142  e->s = -1e10;
143  e->t = -10.0;
144  e->trace = AP->vm_trc;
145 
146  e++;
147  e->c.m = MX_VAL;
148  e->c.mx_val.name = dupstr("VM_MN");
149  e->c.mx_val.p_slope = false;
150  // initialize event values
151  e->s = 1e10;
152  e->t = -10.0;
153  e->trace = AP->vm_trc;
154 
155  if(!get_rank())
156  fprintf(stderr,"Parameter sequence written to action potential statistics file.\n");
157 
158  return 0;
159 }
160 
168 void
170 {
171 
172  if(AP->stats!=NULL)
173  fclose(AP->stats);
174 
175  if(AP->rstats!=NULL)
176  fclose(AP->rstats);
177 
178  if((AP->acts.num_events>0) && (AP->acts.events !=NULL)) {
179  for(int i=0;i<AP->acts.num_events;i++)
180  free(AP->acts.events[i].c.x_thresh.name);
181  free(AP->acts.events);
182  }
183 
184  if((AP->repols.num_events>0) && (AP->repols.events !=NULL)) {
185  for(int i=0;i<AP->repols.num_events;i++)
186  free(AP->repols.events[i].c.x_thresh.name);
187  free(AP->repols.events);
188  }
189 
190  if((AP->ampls.num_events>0) && (AP->ampls.events !=NULL)) {
191  for(int i=0;i<AP->ampls.num_events;i++)
192  free(AP->ampls.events[i].c.x_thresh.name);
193  free(AP->ampls.events);
194  }
195 }
196 
197 
206 void
207 shift_vm_trace(double vm, action_potential *AP, timer_manager *tm)
208 {
209  int vm_hist_len = VM_HIST_LEN;
210  for(int i=1;i<vm_hist_len;i++) {
211  AP->vm_trc [i] = AP->vm_trc [i-1];
212  AP->dvm_trc[i] = AP->dvm_trc[i-1];
213  }
214 
215  AP->vm_trc [0] = vm;
216  AP->dvm_trc[0] = (AP->vm_trc[0]-AP->vm_trc[1])/tm->time_step;
217 }
218 
226 void
228 {
229  double vm_mx = AP->ampls.events[VmMX].s;
230  double vm_mn = AP->ampls.events[VmMN].s;
231  double AP_amplitude = vm_mx-vm_mn;
232  double thresh90 = vm_mn+AP_amplitude*(1.-90./100.);
233  double thresh30 = vm_mn+AP_amplitude*(1.-30./100.);
234 
235  log_msg(NULL, 0, 0,"Measuring AP parameters for calibration." );
236 
237  ap_event *e = AP->acts.events;
238  log_msg(NULL, 0, 0,"Changing activation threshold xAPD90 (%.1f->%.1f).",
239  e[xAPD90].c.x_thresh.threshold,thresh90);
240  e[xAPD90].c.x_thresh.threshold = thresh90;
241 
242  e = AP->repols.events;
243  log_msg(NULL, 0, 0,"Changing repolarisation threshold APD90 (%.1f->%.1f).",
244  e[APD90].c.x_thresh.threshold,thresh90);
245  e[APD90].c.x_thresh.threshold = thresh90;
246 
247  log_msg(NULL, 0, 0,"Changing repolarisation threshold APD30 (%.1f->%.1f).",
248  e[APD30].c.x_thresh.threshold,thresh30);
249  e[APD30].c.x_thresh.threshold = thresh30;
250 }
251 
262 bool
263 check_events(double vm, action_potential *AP, timer_manager *tm)
264 {
265  // add new sample to trace
266  shift_vm_trace(vm,AP,tm);
267 
268  bool act_triggered = false;
269  double prevAct = AP->acts.events->t;
270 
271  for(int i=0; i<AP->acts.num_events; i++) {
272  ap_event *e = AP->acts.events+i;
273  if(e->c.m==X_THRESHOLD)
274  check_threshold(e,tm);
275  else if (e->c.m==MX_RATE)
276  check_mx_rate(e,tm);
277  else if (e->c.m==MX_VAL)
278  check_mx_val(e,tm);
279  else
280  fprintf( stderr, "Not implemented yet.\n" );
281  }
282 
283  for(int i=0; i<AP->repols.num_events; i++) {
284  ap_event *e = AP->repols.events+i;
285  if(e->c.m==X_THRESHOLD)
286  check_threshold(e,tm);
287  else if (e->c.m==MX_RATE)
288  check_mx_rate(e,tm);
289  else
290  fprintf( stderr, "Not implemented yet.\n" );
291  }
292 
293  for(int i=0; i<AP->ampls.num_events; i++) {
294  ap_event *e = AP->ampls.events+i;
295  if (e->c.m==MX_VAL)
296  check_mx_val(e,tm);
297  else
298  fprintf( stderr, "Not implemented yet.\n" );
299  }
300 
301  if(AP->acts.events->trg) {
302 
303  // detect false activation during calibration beats
304  if((AP->beat==AP->calBeat) && (tm->time - prevAct < 5.))
305  return false;
306 
307  AP->APD = AP->repols.events->t-prevAct;
308  AP->DI = tm->time - AP->repols.events->t;
309  AP->triang = AP->repols.events[APD90].t-AP->repols.events[APD30].t;
310 
311  // onset of AP
312  AP->state = AP_STATE;
313  act_triggered = true;
314 
315  AP->beat++;
316  if(AP->beat==AP->calBeat) calibrate_thresholds(AP);
317  if(AP->beat>0) print_AP_stats(AP);
318 
319  AP->DIp = AP->DI;
320 
321  bool is_steady = update_steady_state(AP);
322 
323  log_msg( NULL, 0, 0, "Activation event %s triggered at %.3f ms.",
324  AP->acts.events->c.x_thresh.name, tm->time );
325 
326  // reset min/max search for each AP
327  for(int i=0;i<AP->ampls.num_events;i++) {
328  ap_event *e = AP->ampls.events+i;
329  if(e->c.mx_val.p_slope)
330  e->s = -10e10;
331  else
332  e->s = 10e10;
333  e->t = -10.;
334  }
335  }
336 
337  // entering systolic phase
338  if(AP->acts.events->trg) AP->state = AP_STATE;
339 
340  // entering diastolic phase
341  if(AP->repols.events->trg) AP->state = DI_STATE;
342 
343  return act_triggered;
344 }
345 
354 bool
356 {
357  static int init = 0;
358  bool is_stedy = false;
359 
360  #define NUM_STEADY_APs 5
361  #define SS_PERC_ERR 0.01
362  if(!init) {
363  AP->ss.num_APs = NUM_STEADY_APs;
364  AP->ss.num_params = 5;
365  AP->ss.cnt = 0;
366  AP->ss.c = (float *)calloc(AP->ss.num_params,sizeof(float));
367  AP->ss.p = (float *)calloc(AP->ss.num_params,sizeof(float));
368  init = 1;
369  }
370 
371  // first we update current
372  AP->ss.c[0] = AP->APD;
373  AP->ss.c[1] = AP->DI;
374  AP->ss.c[2] = AP->triang;
375  AP->ss.c[3] = AP->ampls.events[0].s; // maximum Vm over AP
376  AP->ss.c[4] = AP->ampls.events[1].s; // minimum Vm over AP
377 
378  int d = 0;
379  for(int i=0;i<AP->ss.num_params;i++) {
380  float r_err = fabs(AP->ss.c[i]-AP->ss.p[i])/AP->ss.c[i]*100.;
381  if(r_err>SS_PERC_ERR) d++;
382  AP->ss.p[i] = AP->ss.c[i];
383  }
384  if(!d) AP->ss.cnt++;
385  else AP->ss.cnt = 0;
386 
387  return check_steady_state(&AP->ss);
388 }
389 
397 bool
399 {
400  if(ss->cnt>=ss->num_APs)
401  return true;
402  else
403  return false;
404 }
405 
414 bool
415 check_threshold(ap_event *e, timer_manager *tm)
416 {
417  e->trg = false;
418  double *s = e->trace;
419  if(e->c.x_thresh.p_slope) {
420  if(s[0]>=e->c.x_thresh.threshold && s[1]<e->c.x_thresh.threshold) {
421  e->s = e->c.x_thresh.threshold;
422  e->t = tm->time + (e->c.x_thresh.threshold-s[1])/(s[0]-s[1])*tm->time_step;
423  e->trg = true;
424  }
425  }
426  else {
427  if(s[0]<=e->c.x_thresh.threshold && s[1]>e->c.x_thresh.threshold) {
428  e->s = e->c.x_thresh.threshold;
429  e->t = tm->time - (e->c.x_thresh.threshold-s[1])/(s[0]-s[1])*tm->time_step;
430  e->trg = true;
431  }
432  }
433  return e->trg;
434 }
435 
444 bool
445 check_mx_rate(ap_event *e, timer_manager *tm)
446 {
447  e->trg = false;
448  double *s = e->trace;
449  if(e->c.mx_rate.p_slope) {
450  if(s[0]>=e->s) {
451  e->s = s[0];
452  e->t = tm->time - tm->time_step*0.5;
453  e->trg = true;
454  }
455  }
456  else {
457  if(s[0]<e->s) {
458  e->s = s[0];
459  e->t = tm->time - tm->time_step*0.5;
460  e->trg = true;
461  }
462  }
463  return e->trg;
464 }
465 
474 bool
475 check_mx_val(ap_event *e, timer_manager* tm)
476 {
477  e->trg = false;
478  double *s = e->trace;
479  if(e->c.mx_val.p_slope) {
480  if(e->s<s[0]) { e->s = s[0]; e->t = tm->time; }
481  }
482  else
483  if(e->s>s[0]) { e->s = s[0]; e->t = tm->time; }
484 
485  return e->trg;
486 }
487 
497 void
499 {
500  if(get_rank()) return;
501 
502  double vm_mx = AP->ampls.events[0].s;
503  double vm_mn = AP->ampls.events[1].s;
504  double act_tm = AP->acts.events[0].t;
505  bool is_steady = AP->ss.cnt>=AP->ss.num_APs;
506  fprintf(AP->stats,"%d %1s %1d %.2f %.2f %.2f %.2f %.2f %.2f %.3f\n",
507  AP->beat,AP->pmat?"P":"*",is_steady,AP->APD,AP->DI,AP->DIp,
508  AP->triang,vm_mx, vm_mn, act_tm);
509 
510  // running restitution protocol?
511  if(AP->rstats!=NULL && AP->pmat)
512  fprintf(AP->rstats,"%d %1s %1d %.2f %.2f %.2f %.2f %.2f %.2f %.3f\n",
513  AP->beat,AP->pmat?"P":"*",is_steady,AP->APD,AP->DI,AP->DIp,
514  AP->triang,vm_mx, vm_mn, act_tm);
515 }
516 
525 void
527 {
528  if(get_rank()) return;
529  fprintf(outbuf, "# Beat Prematurity(P||*) steady_state ");
530  fprintf(outbuf, "APD(n) DI(n) DI(n-1) Triangulation VmMax VmMin tAct\n" );
531 }
532 
533 } // namespace limpet
void calibrate_thresholds(action_potential *AP)
Definition: ap_analyzer.cc:227
int calBeat
beat number used to calibrate thresholds
Definition: ap_analyzer.h:90
bool check_steady_state(steady_state_ap *ss)
Definition: ap_analyzer.cc:398
bool pmat
indiacte premature AP
Definition: ap_analyzer.h:91
int cnt
count number of beats in steady state
Definition: ap_analyzer.h:80
ap_events ampls
amplitude events (mx or mn found in a signal)
Definition: ap_analyzer.h:94
enum limpet::_act_events act_events
char * dupstr(const char *old_str)
Definition: basics.cc:44
int beat
beat counter
Definition: ap_analyzer.h:89
event_mx_rate mx_rate
Definition: ap_analyzer.h:60
#define NUM_STEADY_APs
float * p
array of parameters at previous AP
Definition: ap_analyzer.h:82
event_mx_val mx_val
Definition: ap_analyzer.h:61
FILE * rstats
output restitution statistics only
Definition: ap_analyzer.h:104
float triang
AP triangulation, APD90-APD30.
Definition: ap_analyzer.h:97
enum limpet::_repol_events repol_events
void print_AP_stats_header(action_potential *AP, FILE *outbuf)
Definition: ap_analyzer.cc:526
double dvm_trc[VM_HIST_LEN]
Definition: ap_analyzer.h:101
bool check_events(double vm, action_potential *AP, timer_manager *tm)
Definition: ap_analyzer.cc:263
double t
time when event was observed
Definition: ap_analyzer.h:66
double * trace
short trace of signal to analyze
Definition: ap_analyzer.h:67
event_x_threshold x_thresh
Definition: ap_analyzer.h:59
bool check_threshold(ap_event *e, timer_manager *tm)
Definition: ap_analyzer.cc:415
double s
value of signal when event occurred
Definition: ap_analyzer.h:65
float * c
array of parameters at current AP
Definition: ap_analyzer.h:81
float APD
standard APD as measured by APD90
Definition: ap_analyzer.h:95
#define SS_PERC_ERR
int initialize_AP_analysis(action_potential *AP)
Definition: ap_analyzer.cc:71
void print_AP_stats(action_potential *AP)
Definition: ap_analyzer.cc:498
#define CALIBRATION_BEAT
Definition: ap_analyzer.h:29
#define def_APD30_thresh
Definition: ap_analyzer.h:33
int num_APs
minimum number of APs for steady state
Definition: ap_analyzer.h:78
void cleanup_AP_analysis(action_potential *AP)
Definition: ap_analyzer.cc:169
double vm_trc[VM_HIST_LEN]
Definition: ap_analyzer.h:100
enum limpet::_ampl_events ampl_events
#define BOGUS
ap_events repols
repolarization events
Definition: ap_analyzer.h:93
bool update_steady_state(action_potential *AP)
Definition: ap_analyzer.cc:355
ap_events acts
activation events
Definition: ap_analyzer.h:92
FILE * stats
output statistics for each AP
Definition: ap_analyzer.h:103
ap_states state
current state of AP, ongoing AP or diastole
Definition: ap_analyzer.h:99
bool p_slope
toggle max min search
Definition: ap_analyzer.h:54
#define def_APD90_thresh
Definition: ap_analyzer.h:32
steady_state_ap ss
steady state detection
Definition: ap_analyzer.h:102
int num_params
number of parameters we use to decide steady state
Definition: ap_analyzer.h:79
void shift_vm_trace(double vm, action_potential *AP, timer_manager *tm)
Definition: ap_analyzer.cc:207
ap_event * events
Definition: ap_analyzer.h:74
void log_msg(FILE_SPEC out, int level, unsigned char flag, const char *fmt,...)
Definition: basics.cc:72
#define VM_HIST_LEN
Definition: ap_analyzer.h:87
float DIp
DI of the previous AP.
Definition: ap_analyzer.h:98
char * name
name of event
Definition: ap_analyzer.h:53
bool check_mx_rate(ap_event *e, timer_manager *tm)
Definition: ap_analyzer.cc:445
bool check_mx_val(ap_event *e, timer_manager *tm)
Definition: ap_analyzer.cc:475
Basic utility structs and functions, mostly IO related.
float DI
BCL-APD90.
Definition: ap_analyzer.h:96
int get_rank(MPI_Comm comm=PETSC_COMM_WORLD)
Definition: basics.h:276
event_class c
Definition: ap_analyzer.h:68