ap_analyzer.cc

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is licensed under the openCARP Academic Public License (APL)
// v1.0: You can use and redistribute it and/or modify it in non-commercial
// academic environments under the terms of APL as published by the openCARP
// project v1.0, or (at your option) any later version. Commercial use requires
// a commercial license (info@opencarp.org).
//
// This program is distributed without any warranty; see the openCARP APL for
// more details.
//
// You should have received a copy of the openCARP APL along with this program
// and can find it online: http://www.opencarp.org/license
// ----------------------------------------------------------------------------

#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "ap_analyzer.h"
#include "basics.h"

namespace limpet {

using ::opencarp::dupstr;
using ::opencarp::get_rank;
using ::opencarp::log_msg;
using ::opencarp::timer_manager;

void shift_vm_trace(double vm, action_potential *AP, timer_manager *tm);
bool check_threshold(ap_event *e, timer_manager *tm);
bool check_mx_rate(ap_event *e,  timer_manager *tm);
bool check_mx_val (ap_event *e, timer_manager *tm);
bool update_steady_state(action_potential *AP);
void print_AP_stats(action_potential *AP);


typedef enum _act_events   {xAPD90,VmDotMx,nActEvents} act_events;
typedef enum _repol_events {APD90,APD30,nRepolEvents} repol_events;
typedef enum _ampl_events  {VmMX,VmMN,nAmplEvents} ampl_events;


int
initialize_AP_analysis(action_potential *AP)
{
  memset(AP,0,sizeof(action_potential));
  AP->beat             = -1;
  // define which beat to use for calibration
  AP->calBeat          = CALIBRATION_BEAT;
  AP->pmat             = false;
  int num_act_events   = nActEvents;
  int num_repol_events = nRepolEvents;
  int num_ampl_events  = nAmplEvents;

#define BOGUS 9999
  // analyze action potentials
  // find activation marker
  AP->acts.num_events     = num_act_events,
  AP->acts.events         = (ap_event *)calloc(num_act_events,sizeof(ap_event));

  // set event type
  ap_event *e             = AP->acts.events;
  e->c.m                  = X_THRESHOLD;
  e->c.x_thresh.name      = dupstr("LAT_xAPD90");
  e->c.x_thresh.threshold = def_APD90_thresh;
  e->c.x_thresh.p_slope   = true;
  // initialize event values
  e->s                    = BOGUS;
  e->t                    = -10.0;
  e->trace                = AP->vm_trc;  // either vm or dvm trace

  e++;
  e->c.m                  = MX_VAL;
  e->c.x_thresh.name      = dupstr("LAT_dVdt");
  e->c.x_thresh.p_slope   = true;
  // initialize event values
  e->s                    = -1e10;
  e->t                    = -10.;
  e->trace                = AP->dvm_trc;  // find max of dvm trace


  // find repolarisation marker
  AP->repols.num_events   = num_repol_events,
  AP->repols.events       = (ap_event *)calloc(num_repol_events,sizeof(ap_event));
  // set event type
  e                       = AP->repols.events;
  e->c.m                  = X_THRESHOLD;
  e->c.x_thresh.name      = dupstr("REPOL_90");
  e->c.x_thresh.threshold = def_APD90_thresh;
  e->c.x_thresh.p_slope   = false;
  // initialize event values
  e->s                    = BOGUS;
  e->t                    = -10.0;
  e->trace                = AP->vm_trc;

  e++;
  e->c.m                  = X_THRESHOLD;
  e->c.x_thresh.name      = dupstr("REPOL_30");
  e->c.x_thresh.threshold = def_APD30_thresh;
  e->c.x_thresh.p_slope   = false;
  // initialize event values
  e->s                    = BOGUS;
  e->t                    = -10.0;
  e->trace                = AP->vm_trc;

  // find extremal signal values
  AP->ampls.num_events    = num_ampl_events,
  AP->ampls.events        = (ap_event *)calloc(num_ampl_events,sizeof(ap_event));
  // set event type
  e                       = AP->ampls.events;
  e->c.m                  = MX_VAL;
  e->c.mx_val.name        = dupstr("VM_MX");
  e->c.mx_val.p_slope     = true;
  // initialize event values
  e->s                    = -1e10;
  e->t                    = -10.0;
  e->trace                = AP->vm_trc;

  e++;
  e->c.m                  = MX_VAL;
  e->c.mx_val.name        = dupstr("VM_MN");
  e->c.mx_val.p_slope     = false;
  // initialize event values
  e->s                    =  1e10;
  e->t                    = -10.0;
  e->trace                = AP->vm_trc;

  if(!get_rank()) 
    fprintf(stderr,"Parameter sequence written to action potential statistics file.\n");
  
  return 0;
}

void
cleanup_AP_analysis(action_potential *AP)
{

  if(AP->stats!=NULL)
    fclose(AP->stats);

  if(AP->rstats!=NULL)
    fclose(AP->rstats);

  if((AP->acts.num_events>0) && (AP->acts.events !=NULL))  {
    for(int i=0;i<AP->acts.num_events;i++)
      free(AP->acts.events[i].c.x_thresh.name);
    free(AP->acts.events);
  }

  if((AP->repols.num_events>0) && (AP->repols.events !=NULL))  {
    for(int i=0;i<AP->repols.num_events;i++)
      free(AP->repols.events[i].c.x_thresh.name);
    free(AP->repols.events);
  }

  if((AP->ampls.num_events>0) && (AP->ampls.events !=NULL))  {
    for(int i=0;i<AP->ampls.num_events;i++)
      free(AP->ampls.events[i].c.x_thresh.name);
    free(AP->ampls.events);
  }
}


void
shift_vm_trace(double vm, action_potential *AP, timer_manager *tm)
{
  int vm_hist_len = VM_HIST_LEN;
  for(int i=1;i<vm_hist_len;i++)  {
    AP->vm_trc [i] = AP->vm_trc [i-1];
    AP->dvm_trc[i] = AP->dvm_trc[i-1];
  }

  AP->vm_trc [0] = vm;
  AP->dvm_trc[0] = (AP->vm_trc[0]-AP->vm_trc[1])/tm->time_step;
}

void
calibrate_thresholds(action_potential *AP)
{
  double vm_mx = AP->ampls.events[VmMX].s;
  double vm_mn = AP->ampls.events[VmMN].s;
  double AP_amplitude = vm_mx-vm_mn;
  double thresh90 = vm_mn+AP_amplitude*(1.-90./100.);
  double thresh30 = vm_mn+AP_amplitude*(1.-30./100.);

  log_msg(NULL, 0, 0,"Measuring AP parameters for calibration." );

  ap_event *e = AP->acts.events;
  log_msg(NULL, 0, 0,"Changing activation threshold xAPD90 (%.1f->%.1f).",
                        e[xAPD90].c.x_thresh.threshold,thresh90);
  e[xAPD90].c.x_thresh.threshold = thresh90;

  e = AP->repols.events;
  log_msg(NULL, 0, 0,"Changing repolarisation threshold APD90 (%.1f->%.1f).",
                        e[APD90].c.x_thresh.threshold,thresh90);
  e[APD90].c.x_thresh.threshold = thresh90;

  log_msg(NULL, 0, 0,"Changing repolarisation threshold APD30 (%.1f->%.1f).",
                        e[APD30].c.x_thresh.threshold,thresh30);
  e[APD30].c.x_thresh.threshold = thresh30;
}

bool
check_events(double vm, action_potential *AP, timer_manager *tm)
{
  // add new sample to trace
  shift_vm_trace(vm,AP,tm);

  bool   act_triggered = false;
  double prevAct       = AP->acts.events->t;

  for(int i=0; i<AP->acts.num_events; i++)  {
    ap_event *e = AP->acts.events+i;
    if(e->c.m==X_THRESHOLD)
      check_threshold(e,tm);
    else if (e->c.m==MX_RATE)
      check_mx_rate(e,tm);
    else if (e->c.m==MX_VAL)
      check_mx_val(e,tm);
    else
      fprintf( stderr, "Not implemented yet.\n" );
  }

  for(int i=0; i<AP->repols.num_events; i++)  {
    ap_event *e = AP->repols.events+i;
    if(e->c.m==X_THRESHOLD)
      check_threshold(e,tm);
    else if (e->c.m==MX_RATE)
      check_mx_rate(e,tm);
    else
      fprintf( stderr, "Not implemented yet.\n" );
  }

  for(int i=0; i<AP->ampls.num_events; i++)  {
    ap_event *e = AP->ampls.events+i;
    if (e->c.m==MX_VAL)
      check_mx_val(e,tm);
    else
      fprintf( stderr, "Not implemented yet.\n" );
  }

  if(AP->acts.events->trg)  {

    // detect false activation during calibration beats
    if((AP->beat==AP->calBeat) && (tm->time - prevAct < 5.))
      return false;

    AP->APD    = AP->repols.events->t-prevAct;
    AP->DI     = tm->time - AP->repols.events->t;
    AP->triang = AP->repols.events[APD90].t-AP->repols.events[APD30].t;

    // onset of AP
    AP->state     = AP_STATE;
    act_triggered = true;

    AP->beat++;
    if(AP->beat==AP->calBeat) calibrate_thresholds(AP);
    if(AP->beat>0)  print_AP_stats(AP);

    AP->DIp = AP->DI;

    bool is_steady = update_steady_state(AP);

    log_msg( NULL, 0, 0, "Activation event %s triggered at %.3f ms.",
             AP->acts.events->c.x_thresh.name, tm->time );

    // reset min/max search for each AP
    for(int i=0;i<AP->ampls.num_events;i++)  {
      ap_event *e = AP->ampls.events+i;
      if(e->c.mx_val.p_slope)
        e->s = -10e10;
      else
        e->s = 10e10;
      e->t = -10.;
    }
  }

  // entering systolic phase
  if(AP->acts.events->trg)   AP->state = AP_STATE;

  // entering diastolic phase
  if(AP->repols.events->trg) AP->state = DI_STATE;

  return act_triggered;
}

bool
update_steady_state(action_potential *AP)
{
  static int init     = 0;
  bool       is_stedy = false;

  #define NUM_STEADY_APs     5
  #define SS_PERC_ERR        0.01
  if(!init)  {
    AP->ss.num_APs    = NUM_STEADY_APs;
    AP->ss.num_params = 5;
    AP->ss.cnt        = 0;
    AP->ss.c          = (float *)calloc(AP->ss.num_params,sizeof(float));
    AP->ss.p          = (float *)calloc(AP->ss.num_params,sizeof(float));
    init              = 1;
  }

  // first we update current
  AP->ss.c[0] = AP->APD;
  AP->ss.c[1] = AP->DI;
  AP->ss.c[2] = AP->triang;
  AP->ss.c[3] = AP->ampls.events[0].s; // maximum Vm over AP
  AP->ss.c[4] = AP->ampls.events[1].s; // minimum Vm over AP

  int d = 0;
  for(int i=0;i<AP->ss.num_params;i++)  {
    float r_err = fabs(AP->ss.c[i]-AP->ss.p[i])/AP->ss.c[i]*100.;
    if(r_err>SS_PERC_ERR) d++;
    AP->ss.p[i] = AP->ss.c[i];
  }
  if(!d)  AP->ss.cnt++;
  else AP->ss.cnt = 0;

  return check_steady_state(&AP->ss);
}

bool
check_steady_state(steady_state_ap *ss)
{
  if(ss->cnt>=ss->num_APs)
    return true;
  else
    return false;
}

bool
check_threshold(ap_event *e, timer_manager *tm)
{
  e->trg = false;
  double *s = e->trace;
  if(e->c.x_thresh.p_slope)  {
    if(s[0]>=e->c.x_thresh.threshold && s[1]<e->c.x_thresh.threshold)  {
      e->s   = e->c.x_thresh.threshold;
      e->t   = tm->time + (e->c.x_thresh.threshold-s[1])/(s[0]-s[1])*tm->time_step;
      e->trg = true;
    }
  }
  else  {
    if(s[0]<=e->c.x_thresh.threshold && s[1]>e->c.x_thresh.threshold)  {
      e->s   = e->c.x_thresh.threshold;
      e->t   = tm->time - (e->c.x_thresh.threshold-s[1])/(s[0]-s[1])*tm->time_step;
      e->trg = true;
    }
  }
  return e->trg;
}

bool
check_mx_rate(ap_event *e, timer_manager *tm)
{
  e->trg = false;
  double *s = e->trace;
  if(e->c.mx_rate.p_slope)  {
    if(s[0]>=e->s)  {
      e->s   = s[0];
      e->t   = tm->time - tm->time_step*0.5;
      e->trg = true;
    }
  }
  else  {
    if(s[0]<e->s)  {
      e->s   = s[0];
      e->t   = tm->time - tm->time_step*0.5;
      e->trg = true;
    }
  }
  return e->trg;
}

bool
check_mx_val(ap_event *e, timer_manager* tm)
{
  e->trg = false;
  double *s = e->trace;
  if(e->c.mx_val.p_slope)  {
    if(e->s<s[0]) { e->s = s[0]; e->t = tm->time; }
  }
  else
    if(e->s>s[0]) { e->s = s[0]; e->t = tm->time; }

  return e->trg;
}

void
print_AP_stats(action_potential *AP)
{
  if(get_rank())  return;

  double vm_mx     = AP->ampls.events[0].s;
  double vm_mn     = AP->ampls.events[1].s;
  double act_tm    = AP->acts.events[0].t;
  bool   is_steady = AP->ss.cnt>=AP->ss.num_APs;
  fprintf(AP->stats,"%d %1s %1d %.2f %.2f %.2f %.2f %.2f %.2f %.3f\n",
          AP->beat,AP->pmat?"P":"*",is_steady,AP->APD,AP->DI,AP->DIp,
          AP->triang,vm_mx, vm_mn, act_tm);

  // running restitution protocol?
  if(AP->rstats!=NULL && AP->pmat)
    fprintf(AP->rstats,"%d %1s %1d %.2f %.2f %.2f %.2f %.2f %.2f %.3f\n",
            AP->beat,AP->pmat?"P":"*",is_steady,AP->APD,AP->DI,AP->DIp,
            AP->triang,vm_mx, vm_mn, act_tm);
}

void
print_AP_stats_header(action_potential *AP, FILE *outbuf)
{
  if(get_rank()) return;
  fprintf(outbuf, "# Beat Prematurity(P||*) steady_state ");
  fprintf(outbuf, "APD(n) DI(n) DI(n-1) Triangulation VmMax VmMin tAct\n" );
}

}  // namespace limpet