sim_av.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpcift.h"
#include "tpctac.h"
#include "tpcpar.h"
#include "tpctacmod.h"
#include "tpccm.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Simulation of venous BTAC (Cv) from arterial BTAC (Ca), based on",
  "three-tissue compartmental model:",
  " ",
  "  ____    K1   ____    k3   ____    k5   ____     ",
  " | Ca | ----> | C1 | ----> | C2 | ----> | C3 |    ",
  " |____| <---- |____| <---- |____| <---- |____|    ",
  "          k2           k4           k6            ",
  " ",
  "  dC1(t)/dt = K1*Ca(T) - (k2+k3)*C1(T) + k4*C2(T) ",
  "  dC2(t)/dt = k3*C1(T) - (k4+k5)*C2(T) + k6*C3(T) ",
  "  dC3(t)/dt = k5*C2(T) - k6*C3(T)                 ",
  "  Ct(T) = C1(T) + C2(T) + C3(T)                   ",
  "  Cv(T) = Ca(T) - dCt(t)/dt / f ",
  " ",
  "Usage: @P [options] parfile [abtacfile simfile]",
  " ",
  "Options:",
  " -ttac=<Filename>",
  "     Simulated tissue curve is written in specified file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "To create a template parameter file, do not enter names for input and",
  "simulated TACs. Obligatory parameters are f, K1, K1/k2, k3, k3/k4.",
  "Parameters k5 and k5/k6 are optional.",
  "If parameter file does not contain units, then per min and per mL units",
  "are assumed.",
  "For accurate results, BTAC should have very short sampling intervals.",
  "To reduce the model, k3 or k3/k4 can be set to 0, to apply one-tissue or",
  "irreversible two-tissue model, respectively, or k5 or k5/k6 can be set to 0,",
  "to apply two-tissue or irreversible three-tissue model, respectively.",
  " ",
  "See also: sim_3tcm, tacadd, taccalc, tac2svg, convexpf, fit_xexp",
  " ",
  "Keywords: input, simulation, modelling, compartmental model",
  0};
/*****************************************************************************/
 
/*****************************************************************************/
/* Turn on the globbing of the command line, since it is disabled by default in
   mingw-w64 (_dowildcard=0); in MinGW32 define _CRT_glob instead, if necessary;
   In Unix&Linux wildcard command line processing is enabled by default. */
/*
#undef _CRT_glob
#define _CRT_glob -1
*/
int _dowildcard = -1;
/*****************************************************************************/
 
/*****************************************************************************/

int simC3vb(
  double *t,
  double *cab,
  const int nr,
  double f,
  double k1,
  double k2,
  double k3,
  double k4,
  double k5,
  double k6,
  double *cvb,
  double *ct
) {
  int i;
  double b, c, d, w, z, dt2;
  double cai, ca_last, t_last;
  double ct1, ct1_last, ct2, ct2_last, ct3, ct3_last;
  double ct1i, ct1i_last, ct2i, ct2i_last, ct3i, ct3i_last;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(cvb==NULL) return 2;
 
  /* Check actual parameter number */
  if(!(f>=0.0) || !(k1>=0.0) || k1>f) return 3;
  if(k3<=0.0) {k3=0.0; if(k2<=0.0) k2=0.0;}
  else if(k5<=0.0) {k5=0.0; if(k4<=0.0) k4=0.0;}
  else {if(k6<=0.0) k6=0.0;}
 
  /* Calculate curves */
  t_last=0.0; if(t[0]<t_last) t_last=t[0]; 
  cai=ca_last=0.0;
  ct1_last=ct2_last=ct3_last=ct1i_last=ct2i_last=ct3i_last=0.0;
  ct1=ct2=ct3=ct1i=ct2i=ct3i=0.0;
  for(i=0; i<nr; i++) {
    /* delta time / 2 */
    dt2=0.5*(t[i]-t_last);
    /* calculate values */
    if(dt2<0.0) {
      return 5;
    } else if(dt2>0.0) {
      /* arterial integral */
      cai+=(cab[i]+ca_last)*dt2;
      /* partial results */
      b=ct1i_last+dt2*ct1_last;
      c=ct2i_last+dt2*ct2_last;
      d=ct3i_last+dt2*ct3_last;
      w=k4 + k5 - (k5*k6*dt2)/(1.0+k6*dt2);
      z=1.0+w*dt2;
      /* 1st tissue compartment and its integral */
      ct1 = (
          + k1*z*cai + (k3*k4*dt2 - (k2+k3)*z)*b
          + k4*c + k4*k6*dt2*d/(1.0+k6*dt2)
        ) / ( z*(1.0 + dt2*(k2+k3)) - k3*k4*dt2*dt2 );
      ct1i = ct1i_last + dt2*(ct1_last+ct1);
      /* 2nd tissue compartment and its integral */
      ct2 = (k3*ct1i - w*c + k6*d/(1.0+k6*dt2)) / z;
      ct2i = ct2i_last + dt2*(ct2_last+ct2);
      /* 3rd tissue compartment and its integral */
      ct3 = (k5*ct2i - k6*d) / (1.0 + k6*dt2);
      ct3i = ct3i_last + dt2*(ct3_last+ct3);
    }
    if(f>0.0) {
      double dct = k1*cab[i] - k2*ct1; 
      cvb[i] = cab[i] - dct/f;
    } else
      cvb[i]=cab[i];
 
    /* copy values to argument arrays; set very small values to zero */
    if(ct!=NULL) {ct[i]=ct1+ct2+ct3; if(fabs(ct[i])<1.0e-12) ct[i]=0.0;}
 
    /* prepare to the next loop */
    t_last=t[i]; ca_last=cab[i];
    ct1_last=ct1; ct1i_last=ct1i;
    ct2_last=ct2; ct2i_last=ct2i;
    ct3_last=ct3; ct3i_last=ct3i;
  }
 
  return 0;
}

/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int    ai, help=0, version=0, verbose=1;
  char   blofile[FILENAME_MAX], simfile[FILENAME_MAX], parfile[FILENAME_MAX];
  char   tisfile[FILENAME_MAX];
  unsigned int parNr=0;
  PAR    par;
  int    ret;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  blofile[0]=simfile[0]=parfile[0]=tisfile[0]=(char)0;
  parInit(&par);
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(!*cptr) continue;
    if(strncasecmp(cptr, "TTAC=", 5)==0 && strlen(cptr)>5) {
      strlcpy(tisfile, cptr+5, FILENAME_MAX); continue;
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break; // tac name argument may start with '-'
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-1;
  
  /* Print help or version? */
  if(help==2) {tpcHtmlUsage(argv[0], info, ""); return(0);}
  if(help) {tpcPrintUsage(argv[0], info, stdout); return(0);}
  if(version) {tpcPrintBuild(argv[0], stdout); return(0);}
 
  /* Process other arguments, starting from the first non-option */
  if(ai<argc) strlcpy(parfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(blofile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(simfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!parfile[0]) {
    fprintf(stderr, "Error: missing parameter file; use option --help\n");
    return(1);
  }
  if(blofile[0] && !simfile[0]) {
    fprintf(stderr, "Error: missing filename.\n");
    return(1);
  }
  if(tisfile[0] && (strcasecmp(tisfile, blofile)==0 || strcasecmp(tisfile, simfile)==0)) {
    fprintf(stderr, "Error: invalid file name for simulated TTAC.\n");
    return(1);
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("parfile := %s\n", parfile);
    printf("abtacfile := %s\n", blofile);
    printf("simfile := %s\n", simfile);
    if(!tisfile[0]) printf("ttacfile := %s\n", tisfile);
  }
  
  /*
   *  Make template parameter file, if no other file names were given, and then quit
   */
  if(!blofile[0]) {
    parNr=7;
    ret=parAllocate(&par, parNr, 3);
    if(ret) {
      fprintf(stderr, "Error: cannot allocate memory for parameters.\n");
      parFree(&par); return(1);
    }
    par.parNr=parNr; par.tacNr=3;
    strcpy(par.n[0].name, "f"); par.n[0].unit=UNIT_ML_PER_ML_MIN;
    strcpy(par.n[1].name, "K1"); par.n[1].unit=UNIT_ML_PER_ML_MIN;
    strcpy(par.n[2].name, "K1/k2"); par.n[2].unit=UNIT_ML_PER_ML;
    strcpy(par.n[3].name, "k3"); par.n[3].unit=UNIT_PER_MIN;
    strcpy(par.n[4].name, "k3/k4"); par.n[4].unit=UNIT_UNITLESS;
    strcpy(par.n[5].name, "k5"); par.n[5].unit=UNIT_PER_MIN;
    strcpy(par.n[6].name, "k5/k6"); par.n[6].unit=UNIT_UNITLESS;
    for(int i=0; i<par.tacNr; i++) {
      sprintf(par.r[i].name, "btac%d", 1+i);
      par.r[i].model=0;
      par.r[i].p[0]=0.2;
      par.r[i].p[1]=0.1;
      par.r[i].p[2]=1.0;
      par.r[i].p[3]=0.05;
      par.r[i].p[4]=1.0;
      par.r[i].p[5]=0.0;
      par.r[i].p[6]=0.0;
    }
    int i=1;
    par.r[i].p[1]=0.1;
    par.r[i].p[2]=1.0;
    par.r[i].p[3]=0.1;
    par.r[i].p[4]=2.0;
    i=2;
    par.r[i].p[1]=0.1;
    par.r[i].p[2]=1.0;
    par.r[i].p[3]=0.1;
    par.r[i].p[4]=0.5;
    if(verbose>1) printf("writing %s\n", parfile);
    FILE *fp; fp=fopen(parfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing (%s)\n", parfile);
      parFree(&par); return(11);
    }
    ret=parWrite(&par, fp, PAR_FORMAT_TSV_UK, 1, &status);
    fclose(fp); parFree(&par); 
    if(ret!=TPCERROR_OK) {fprintf(stderr, "Error: cannot write %s\n", parfile); return(11);}
    return(0);
  }
 
  
  /*
   *  Read parameter file to get parameters for the simulation, and 
   *  set model in case it was given with command-line option.
   */
  if(verbose>1) fprintf(stdout, "reading %s\n", parfile);
  if(parRead(&par, parfile, &status)) {
    fprintf(stderr, "Error: %s (%s)\n", errorMsg(status.error), parfile);
    parFree(&par);
    return(3);
  }
  if(verbose>5) parWrite(&par, stdout, PAR_FORMAT_TSV_UK, 1, NULL);
  /* Check the validity of model and its parameters */
  for(int i=0; i<par.tacNr; i++) {
    /* check that we got at least five parameters */
    if((unsigned int)par.parNr<5) {
      fprintf(stderr, "Error: invalid parameters for selected model.\n");
      parFree(&par); return(3);
    }
    /* check that we can find the obligatory model parameters */
    ret=0;
    if(parFindParameter(&par, "f")<0) ret++;
    if(parFindParameter(&par, "K1")<0) ret++;
    if(parFindParameter(&par, "K1/k2")<0) ret++;
    if(parFindParameter(&par, "k3")<0) ret++;
    if(parFindParameter(&par, "k3/k4")<0) ret++;
    if(ret) {
      fprintf(stderr, "Error: required parameters not available.\n");
      parFree(&par); return(3);
    }
  }
 
  /*
   *  Read input TACs
   */
  if(verbose>1) fprintf(stdout, "reading input TAC\n");
  TAC input; tacInit(&input);
  ret=tacReadModelingInput(blofile, NULL, NULL, &input, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s (input files)\n", errorMsg(status.error));
    tacFree(&input); parFree(&par); return(4);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(input.format));
    printf("tacNr := %d\n", input.tacNr);
    printf("sampleNr := %d\n", input.sampleNr);
    printf("xunit := %s\n", unitName(input.tunit));
    printf("yunit := %s\n", unitName(input.cunit));
  }
  if(verbose>10) tacWrite(&input, stdout, TAC_FORMAT_PMOD, 1, NULL);
  if(input.sampleNr<3) {
    fprintf(stderr, "Error: too few samples in input TAC.\n");
    tacFree(&input); parFree(&par); return(4);
  }
  if(input.sampleNr<10) {
    fprintf(stderr, "Warning: too few samples for reliable simulation.\n");
  }
 
 
  /*
   *  Check parameter and data units
   */
  {
    int i, punit;
    /* Perfusion */
    i=parFindParameter(&par, "f");
    if(i>=0) {
      punit=par.n[i].unit;
      if(punit==UNIT_ML_PER_DL_MIN) { 
        if(verbose>1) printf("Note: converting f from per dL to per mL.\n");
        par.n[i].unit=UNIT_ML_PER_ML_MIN;
        for(int j=0; j<par.tacNr; j++) par.r[j].p[i]*=0.01;
      }
      if(input.tunit==UNIT_MIN && (punit==UNIT_PER_SEC || punit==UNIT_ML_PER_ML_SEC)) {
        if(verbose>1) printf("Note: converting f from per sec to per min.\n");
        if(punit==UNIT_PER_SEC) par.n[i].unit=UNIT_PER_MIN;
        else par.n[i].unit=UNIT_ML_PER_ML_MIN;
        for(int j=0; j<par.tacNr; j++) par.r[j].p[i]*=60.;
      } else if(input.tunit==UNIT_SEC && (punit==UNIT_PER_MIN || punit==UNIT_ML_PER_ML_MIN)) {
        if(verbose>1) printf("Note: converting f from per min to per sec.\n");
        if(punit==UNIT_PER_MIN) par.n[i].unit=UNIT_PER_SEC;
        else par.n[i].unit=UNIT_ML_PER_ML_SEC;
        for(int j=0; j<par.tacNr; j++) par.r[j].p[i]/=60.;
      }
    }
    /* Rate constants */
    for(int ri=1; ri<=6; ri++) {
      char rcname[3]; sprintf(rcname, "k%d", ri);
      if((i=parFindParameter(&par, rcname))>=0) {
        punit=par.n[i].unit;
        if(input.tunit==UNIT_MIN && (punit==UNIT_PER_SEC || punit==UNIT_ML_PER_ML_SEC)) {
          if(verbose>1) printf("Note: converting %s from per sec to per min.\n", rcname);
          if(punit==UNIT_PER_SEC) par.n[i].unit=UNIT_PER_MIN;
          else par.n[i].unit=UNIT_ML_PER_ML_MIN;
          for(int j=0; j<par.tacNr; j++) par.r[j].p[i]*=60.;
        } else if(input.tunit==UNIT_SEC && (punit==UNIT_PER_MIN || punit==UNIT_ML_PER_ML_MIN)) {
          if(verbose>1) printf("Note: converting %s from per min to per sec.\n", rcname);
          if(punit==UNIT_PER_MIN) par.n[i].unit=UNIT_PER_SEC;
          else par.n[i].unit=UNIT_ML_PER_ML_SEC;
          for(int j=0; j<par.tacNr; j++) par.r[j].p[i]/=60.;
        }
      }
    }
  }
 
 
  /*
   *  Allocate space for simulated data
   */
  if(verbose>1) fprintf(stdout, "allocating space for simulated BTACs\n");
  TAC sim; tacInit(&sim);
  ret=tacAllocateWithPAR(&sim, &par, input.sampleNr, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&input); tacFree(&sim); parFree(&par); return(4);
  }
  sim.format=TAC_FORMAT_PMOD;
  sim.weighting=WEIGHTING_OFF;
  sim.cunit=input.cunit;
  sim.tunit=input.tunit;
  sim.format=input.format;
  sim.isframe=0;
  for(int j=0; j<sim.sampleNr; j++) sim.x[j]=input.x[j];
 
  TAC ttac; tacInit(&ttac);
  if(tisfile[0]) {
    if(verbose>1) fprintf(stdout, "allocating space for simulated TTACs\n");
    ret=tacDuplicate(&sim, &ttac);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&sim); tacFree(&input); tacFree(&sim); parFree(&par); return(4);
    }
  }
 
 
  /*
   *  Simulation
   */
  if(verbose>1) printf("simulating\n");
  double *px=input.x;
  double *py=input.c[0].y;
  double K1, k2, k3, k4, k5, k6, Flow;
  int i;
  for(i=0; i<par.tacNr; i++) {
    if(verbose>2) printf("simulating %s\n", sim.c[i].name);
    /* Create pointers for simulated data. */
    double *pcvb=sim.c[i].y;
    double *pct=NULL; if(tisfile[0]) pct=ttac.c[i].y;
    /* Get parameters */
    K1=k2=k3=k4=Flow=nan("");
    Flow=parGetParameter(&par, "f", i);
    if(isnan(Flow)) Flow=parGetParameter(&par, "Flow", i);
    K1=parGetParameter(&par, "K1", i);
    k2=parGetParameter(&par, "k2", i);
    if(isnan(k2)) {
      double r=parGetParameter(&par, "K1/k2", i);
      if(!(r>0.0)) k2=nan(""); else k2=K1/r;
    }
    k3=parGetParameter(&par, "k3", i);
    k4=parGetParameter(&par, "k4", i);
    if(isnan(k4)) {
      double r=parGetParameter(&par, "k3/k4", i);
      if(!(r>0.0)) k4=nan(""); else k4=k3/r;
    }
    k5=parGetParameter(&par, "k5", i);
    k6=parGetParameter(&par, "k6", i);
    if(isnan(k6)) {
      double r=parGetParameter(&par, "k5/k6", i);
      if(!(r>0.0)) k6=nan(""); else k6=k5/r;
    }
    /* Verify and fix parameters */
    if(!(Flow>=0.0) || !(K1>=0.0) || k2<0 || k3<0 || k4<0 || k5<0 || k6<0) {
      fprintf(stderr, "Error: invalid rate constant.\n");
      tacFree(&sim); tacFree(&ttac); tacFree(&input); parFree(&par); return(5);
    }
    if(K1>Flow) {
      fprintf(stderr, "Error: K1 cannot be higher than f.\n");
      tacFree(&sim); tacFree(&ttac); tacFree(&input); parFree(&par); return(5);
    }
    if(Flow==0.0) {K1=k2=k3=k4=k5=k6=0.0;}
    if(K1==0.0) {k2=k3=k4=k5=k6=0.0;}
    if(k3==0.0) {k4=k5=k6=0.0;}
    if(k5==0.0) k6=0.0;
    if(verbose>3) {
      printf("f := %g\n", Flow);
      printf("K1 := %g\n", K1);
      printf("k2 := %g\n", k2);
      if(!isnan(k3)) printf("k3 := %g\n", k3);
      if(!isnan(k4)) printf("k4 := %g\n", k4);
      if(!isnan(k5)) printf("k5 := %g\n", k5);
      if(!isnan(k6)) printf("k6 := %g\n", k6);
    }
    /* Simulate */
    if(isnan(Flow)) Flow=0.0;
    if(isnan(K1)) K1=0.0; 
    if(isnan(k2)) k2=0.0;
    if(isnan(k3)) k3=0.0; 
    if(isnan(k4)) k4=0.0;
    if(isnan(k5)) k5=0.0;
    if(isnan(k6)) k6=0.0;
    ret=simC3vb(px, py, sim.sampleNr, Flow, K1, k2, k3, k4, k5, k6, pcvb, pct);
    if(ret) {
      fprintf(stderr, "Error: invalid data for simulation.\n");
      if(verbose>1) printf("sim_return_code := %d\n", ret);
      tacFree(&sim); tacFree(&ttac); tacFree(&input); parFree(&par); return(6);
    }
  }
 
  /* simulation done */
  tacFree(&input); parFree(&par);
 
 
  /*
   *  Save simulated data 
   */
  if(verbose>1) printf("writing %s\n", simfile);
  FILE *fp; fp=fopen(simfile, "w");
  if(fp==NULL) {
    fprintf(stderr, "Error: cannot open file for writing (%s)\n", simfile);
    tacFree(&sim); tacFree(&ttac); return(11);
  }
  ret=tacWrite(&sim, fp, TAC_FORMAT_PMOD, 1, &status);
  fclose(fp); tacFree(&sim);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
    tacFree(&ttac); return(12);
  }
  if(verbose>=0) printf("%s saved.\n", simfile);
 
  if(tisfile[0]) {
    if(verbose>1) printf("writing %s\n", tisfile);
    FILE *fp; fp=fopen(tisfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing (%s)\n", tisfile);
      tacFree(&ttac); return(13);
    }
    ret=tacWrite(&ttac, fp, TAC_FORMAT_PMOD, 1, &status);
    fclose(fp); tacFree(&ttac);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
      return(14);
    }
    if(verbose>=0) printf("%s saved.\n", tisfile);
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/