sim_wliv.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 "tpcfileutil.h"
#include "tpcift.h"
#include "tpctac.h"
#include "tpcpar.h"
#include "tpctacmod.h"
#include "tpccm.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Simulation of PET tissue time-radioactivity concentration curves (TACs)",
  "from arterial blood (Ca) TACs, based on compartmental model for radiowater",
  "kinetics in the liver. Liver has two inputs, hepatic artery and portal vein.",
  "Portal vein input is simulated by assuming that GI tract can be modelled",
  "with time delay and dispersion of the arterial input function.",
  "Model has two parameters for the perfusion, represented by K1a and K1p for",
  "the arterial and portal blood flow, respectively, and a single parameter p",
  "for the partition coefficient of water (p).",
  "There are two delay times, LdT, which is the delay between the Ca and liver,",
  "and PdT, which is the additional delay in the portal vein blood.",
  "Dispersion the the GI tract is accounted for by rate constant kGI.",
  "The volume fraction of arterial blood (Vb) is assumed to contain arterial",
  "and portal blood in proportion to their contribution to the total blood flow.",
  " ",
  "Usage: @P [options] parfile [bloodfile simfile]",
  " ",
  "Options:",
//  " -TTM",
//  "     Use transit-time model for simulating portal vein curve.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "To create a template parameter file, do not enter names for arterial input",
  "and simulated TACs.",
  "If parameter file does not contain units, then per min and per mL units",
  "are assumed for K1a, K1p, p, kGI, and Vb, and seconds for deley times.",
  "For accurate results, blood TAC should have very short sampling intervals.",
  " ",
  "See also: sim_h2o, fit_wliv, liverpv, simdisp, fit_h2o, tacadd, simframe",
  " ",
  "Keywords: liver, simulation, compartmental model, perfusion, radiowater",
  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 main(int argc, char **argv)
{
  int    ai, help=0, version=0, verbose=1;
  char   blofile[FILENAME_MAX], simfile[FILENAME_MAX], parfile[FILENAME_MAX];
  PAR    par;
  int    TTM=0; // Use transit-time model instead of portal compartment (1) or not (0)
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  blofile[0]=simfile[0]=parfile[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(strcasecmp(cptr, "TTM")==0) {
      TTM=1; 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]) {
    if(!simfile[0]) {fprintf(stderr, "Error: missing file name.\n"); return(1);}
    if(!fileExist(parfile)) {fprintf(stderr, "Error: parameter file not found.\n"); return(1);}
  } else {
    if(fileExist(parfile)) {fprintf(stderr, "Error: parameter template exists.\n"); return(1);}
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("parfile := %s\n", parfile);
    printf("blofile := %s\n", blofile);
    printf("simfile := %s\n", simfile);
    printf("TTM := %d\n", TTM);
    fflush(stdout);
  }
 
 
  /*
   *  Make template parameter file, if no other file names were given, and then quit
   */
  if(!blofile[0]) {
    /* We have already checked that parameter file does not exist */
 
     /* Allocate space for 6 parameters, 4 TTACs */
    if(parAllocate(&par, 7, 3)) {
      fprintf(stderr, "Error: cannot allocate memory for parameters.\n");
      parFree(&par); return(1);
    }
    par.parNr=7; par.tacNr=3;
    strcpy(par.n[0].name, "K1a"); par.n[0].unit=UNIT_ML_PER_ML_MIN;
    strcpy(par.n[1].name, "K1p"); par.n[1].unit=UNIT_ML_PER_ML_MIN;
    strcpy(par.n[2].name, "p");   par.n[2].unit=UNIT_ML_PER_ML;
    strcpy(par.n[3].name, "Vb");  par.n[3].unit=UNIT_ML_PER_ML;
    strcpy(par.n[4].name, "LdT"); par.n[4].unit=UNIT_SEC;
    strcpy(par.n[5].name, "PdT"); par.n[5].unit=UNIT_SEC;
    strcpy(par.n[6].name, "kGI"); par.n[6].unit=UNIT_ML_PER_ML_MIN;
    for(int i=0; i<par.tacNr; i++) {
      sprintf(par.r[i].name, "tac%d", 1+i);
      par.r[i].model=modelCodeIndex("radiowater-liver");
    }
    int ri=0; // first TTAC
    par.r[ri].p[0]=0.5;
    par.r[ri].p[1]=3.5;
    par.r[ri].p[2]=0.8;
    par.r[ri].p[3]=0.02;
    par.r[ri].p[4]=10.0;
    par.r[ri].p[5]=15.0;
    par.r[ri].p[6]=1.0;
    ri++;
    par.r[ri].p[0]=0.4;
    par.r[ri].p[1]=2.6;
    par.r[ri].p[2]=0.6;
    par.r[ri].p[3]=0.02;
    par.r[ri].p[4]=10.0;
    par.r[ri].p[5]=15.0;
    par.r[ri].p[6]=1.0;
    ri++;
    par.r[ri].p[0]=0.8;
    par.r[ri].p[1]=0.0;
    par.r[ri].p[2]=0.8;
    par.r[ri].p[3]=0.02;
    par.r[ri].p[4]=10.0;
    par.r[ri].p[5]=15.0;
    par.r[ri].p[6]=1.0;
 
    if(TTM) {
      /* In transit-time model PdT and kGI are replaced by beta */
      par.parNr=6;
      strcpy(par.n[5].name, "beta"); par.n[5].unit=UNIT_SEC;
      for(ri=0; ri<par.tacNr; ri++) {
        par.r[ri].model=modelCodeIndex("radiowater-liver-TTM");
        par.r[ri].p[5]=60.0*2.17;
      }
    }
 
    if(verbose>1) {printf("writing %s\n", parfile); fflush(stdout);}
 
    FILE *fp; fp=fopen(parfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing (%s)\n", parfile); fflush(stderr);
      parFree(&par); return(11);
    }
    int 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(12);
    }
    if(verbose>=0) {printf("%s saved.\n", parfile); fflush(stdout);}
    return(0);
  }
 
  
  /*
   *  Read parameter file to get parameters for the simulation
   */
  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 */
  unsigned int model=par.r[0].model;
  if(model!=modelCodeIndex("radiowater-liver") && model!=modelCodeIndex("radiowater-liver-TTM")) {
    fprintf(stderr, "Error: invalid model in parameter file.\n");
    parFree(&par); return(3);
  }
  for(int i=0; i<par.tacNr; i++) {
    if(verbose>3) printf("model %s for tac %s\n", modelCode(par.r[i].model), par.r[i].name);
    if(par.r[i].model!=model) {
      fprintf(stderr, "Error: different models in parameter file.\n");
      parFree(&par); return(3);
    }
    /* check that we got at least as many parameters as the model has */
    if(par.parNr<(int)modelParNr(par.r[i].model)) {
      fprintf(stderr, "Error: invalid parameters for selected model.\n");
      parFree(&par); return(3);
    }
  }
  /* check that we can find the obligatory model parameters */
  int i_LdT=-1, i_PdT=-1, i_kGI=-1, i_K1a=-1, i_K1p=-1, i_p=-1, i_Vb=-1, i_beta=-1;
  if(model==modelCodeIndex("radiowater-liver")) {
    int ret=0;
    if((i_K1a=parFindParameter(&par, "K1a"))<0) ret++;
    if((i_K1p=parFindParameter(&par, "K1p"))<0) ret++;
    if((i_p=parFindParameter(&par, "p"))<0) ret++;
    if((i_Vb=parFindParameter(&par, "Vb"))<0) ret++;
    if((i_LdT=parFindParameter(&par, "LdT"))<0) ret++;
    if((i_PdT=parFindParameter(&par, "PdT"))<0) ret++;
    if((i_kGI=parFindParameter(&par, "kGI"))<0) ret++;
    if(ret) {
      fprintf(stderr, "Error: required parameters not available.\n");
      parFree(&par); return(3);
    }
    if(verbose>5) {
      printf("parameter indices:\n");
      printf("  K1a=p[%d]\n", i_K1a);
      printf("  K1p=p[%d]\n", i_K1p);
      printf("  p=p[%d]\n", i_p);
      printf("  Vb=p[%d]\n", i_Vb);
      printf("  LdT=p[%d]\n", i_LdT);
      printf("  PdT=p[%d]\n", i_PdT);
      printf("  kGI=p[%d]\n", i_kGI);
    }
  } else { // Transit-time model
    int ret=0;
    if((i_K1a=parFindParameter(&par, "K1a"))<0) ret++;
    if((i_K1p=parFindParameter(&par, "K1p"))<0) ret++;
    if((i_p=parFindParameter(&par, "p"))<0) ret++;
    if((i_Vb=parFindParameter(&par, "Vb"))<0) ret++;
    if((i_LdT=parFindParameter(&par, "LdT"))<0) ret++;
    if((i_beta=parFindParameter(&par, "beta"))<0) ret++;
    if(ret) {
      fprintf(stderr, "Error: required parameters not available.\n");
      parFree(&par); return(3);
    }
    if(verbose>5) {
      printf("parameter indices:\n");
      printf("  K1a=p[%d]\n", i_K1a);
      printf("  K1p=p[%d]\n", i_K1p);
      printf("  p=p[%d]\n", i_p);
      printf("  Vb=p[%d]\n", i_Vb);
      printf("  LdT=p[%d]\n", i_LdT);
      printf("  beta=p[%d]\n", i_beta);
    }
  }
 
 
  /*
   *  Read input TAC
   */
  if(verbose>1) fprintf(stdout, "reading input TAC\n");
  TAC input; tacInit(&input);
  if(tacReadModelingInput(blofile, NULL, NULL, &input, &status)!=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);
  /* Check for missing sample times */
  if(tacXNaNs(&input)>0) {
    fprintf(stderr, "Error: missing sample times.\n");
    tacFree(&input); parFree(&par); return(4);
  }
  /* Check for missing concentrations */
  if(tacYNaNs(&input, -1)>0) {
    fprintf(stderr, "Error: missing concentrations.\n");
    tacFree(&input); parFree(&par); return(4);
  }
  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"); fflush(stderr);
  }
  if(tacXUnitConvert(&input, UNIT_SEC, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&input); parFree(&par); return(4);
  }
  if(tacSortByTime(&input, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: invalid sample times.\n");
    tacFree(&input); parFree(&par); return(4);
  }
 
 
  /*
   *  Allocate space for simulated data
   */
  if(verbose>1) fprintf(stdout, "allocating space for simulated TACs\n");
  TAC sim; tacInit(&sim);
  if(tacAllocateWithPAR(&sim, &par, input.sampleNr, &status)!=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];
 
 
  /*
   *  Simulation
   */
  if(verbose>1) {printf("simulating...\n"); fflush(stdout);}
  for(int ri=0; ri<sim.tacNr; ri++) {
    if(verbose>2) {printf("  %s\n", sim.c[ri].name); fflush(stdout);}
    double cf;
 
    /* Make delayed arterial input TAC */
    double LdT=par.r[ri].p[i_LdT];
    cf=unitConversionFactor(par.n[i_LdT].unit, UNIT_SEC); if(isfinite(cf)) LdT*=cf;
    if(verbose>4) printf("    LdT := %g s\n", LdT);
    double ca[sim.sampleNr];
    if(fabs(LdT)<1.0E-08) {
      for(int i=0; i<sim.sampleNr; i++) ca[i]=input.c[0].y[i];
    } else {
      double buf[sim.sampleNr]; for(int i=0; i<sim.sampleNr; i++) buf[i]=input.x[i]+LdT;
      if(liInterpolate(buf, input.c[0].y, sim.sampleNr, sim.x, ca, NULL, NULL, sim.sampleNr, 3, 1, 0))
      {
        fprintf(stderr, "Error: cannot simulate arterial delay.\n");
        tacFree(&input); tacFree(&sim); parFree(&par); return(5);
      }
    }
 
    /* Make delayed and dispersed portal vein TAC */
    double cp[sim.sampleNr]; for(int i=0; i<sim.sampleNr; i++) cp[i]=0.0;
    if(model==modelCodeIndex("radiowater-liver")) {
      double PdT=par.r[ri].p[i_PdT];
      cf=unitConversionFactor(par.n[i_PdT].unit, UNIT_SEC); if(isfinite(cf)) PdT*=cf;
      if(verbose>4) printf("    PdT := %g s\n", PdT);
      if(fabs(PdT)<1.0E-08) {
        for(int i=0; i<sim.sampleNr; i++) cp[i]=ca[i]; // including LdT
      } else {
        double buf[sim.sampleNr]; for(int i=0; i<sim.sampleNr; i++) buf[i]=input.x[i]+LdT+PdT;
        if(liInterpolate(buf, input.c[0].y, sim.sampleNr, sim.x, cp, NULL, NULL, sim.sampleNr, 3, 1, 0))
        {
          fprintf(stderr, "Error: cannot simulate portal delay.\n");
          tacFree(&input); tacFree(&sim); parFree(&par); return(5);
        }
      }
      double kGI=par.r[ri].p[i_kGI];
      cf=unitConversionFactor(par.n[i_kGI].unit, UNIT_ML_PER_ML_SEC); 
      if(!isfinite(cf)) cf=unitConversionFactor(par.n[i_kGI].unit, UNIT_PER_SEC); 
      if(isfinite(cf)) kGI*=cf; else kGI/=60.0;
      if(verbose>4) printf("    kGI := %g 1/s\n", kGI);
      double tau=0.0; if(kGI>0.0) tau=1.0/kGI;
      if(verbose>4) printf("    tauGI := %g s\n", tau);
      if(simDispersion(sim.x, cp, sim.sampleNr, tau, 0.0, NULL)) {
        fprintf(stderr, "Error: cannot simulate dispersion.\n");
        tacFree(&input); tacFree(&sim); parFree(&par); return(5);
      }
    } else { // Transit-time model
      double beta=par.r[ri].p[i_beta];
      cf=unitConversionFactor(par.n[i_beta].unit, UNIT_SEC); if(isfinite(cf)) beta*=cf;
      if(verbose>4) printf("    beta := %g s\n", beta);
 
      // To be continued...
 
 
 
    }
 
    /* Simulate pure tissue TAC */
    double K1a=par.r[ri].p[i_K1a];
    cf=unitConversionFactor(par.n[i_K1a].unit, UNIT_ML_PER_ML_SEC); 
    if(isfinite(cf)) K1a*=cf; else K1a/=60.0;
    if(verbose>4) printf("    K1a := %g mL/(s*mL)\n", K1a);
    if(!(K1a>=0.0)) {
      fprintf(stderr, "Error: invalid K1a for '%s'.\n", sim.c[ri].name);
      tacFree(&input); tacFree(&sim); parFree(&par); return(5);
    }
    double K1p=par.r[ri].p[i_K1p];
    cf=unitConversionFactor(par.n[i_K1p].unit, UNIT_ML_PER_ML_SEC); 
    if(isfinite(cf)) K1p*=cf; else K1p/=60.0;
    if(verbose>4) printf("    K1p := %g mL/(s*mL)\n", K1p);
    if(!(K1p>=0.0)) {
      fprintf(stderr, "Error: invalid K1p for '%s'.\n", sim.c[ri].name);
      tacFree(&input); tacFree(&sim); parFree(&par); return(5);
    }
    double p=par.r[ri].p[i_p];
    if(!(p>0.0) && (K1a+K1p)>0.0) {
      fprintf(stderr, "Error: invalid p for '%s'.\n", sim.c[ri].name);
      tacFree(&input); tacFree(&sim); parFree(&par); return(5);
    }
    cf=unitConversionFactor(par.n[i_p].unit, UNIT_ML_PER_ML); if(isfinite(cf)) p*=cf;
    if(verbose>4) printf("    p := %g mL/mL\n", p);
    double k2=0.0; 
    if(p>0.0) {
      k2=(K1a+K1p)/p; if(!(k2>=0.0)) {
        fprintf(stderr, "Error: invalid k2 for '%s'.\n", sim.c[ri].name);
        tacFree(&input); tacFree(&sim); parFree(&par); return(5);
      }
    }
    if(verbose>4) printf("    k2 := %g 1/s\n", k2);
 
    if(simC1DI(sim.x, ca, cp, sim.sampleNr, K1a, K1p, k2, sim.c[ri].y)) {
      fprintf(stderr, "Error: cannot simulate '%s'.\n", sim.c[ri].name);
      tacFree(&input); tacFree(&sim); parFree(&par); return(5);
    }
 
    /* Add the contribution of arterial blood in tissue vasculature */
    double Vb=par.r[ri].p[i_Vb];
    cf=unitConversionFactor(par.n[i_Vb].unit, UNIT_ML_PER_ML); if(isfinite(cf)) Vb*=cf;
    if(verbose>4) printf("    Vb := %g mL/mL\n", Vb);
    if(!(Vb>=0.0) || Vb>1.0) {
      fprintf(stderr, "Error: invalid Vb for '%s'.\n", sim.c[ri].name);
      tacFree(&input); tacFree(&sim); parFree(&par); return(5);
    }
    if(Vb>0.0) {
      for(int i=0; i<sim.sampleNr; i++) sim.c[ri].y[i]*=(1.0-Vb);
      if(K1a+K1p > 0.0)
        for(int i=0; i<sim.sampleNr; i++) sim.c[ri].y[i]+=Vb*(K1a*ca[i]+K1p*cp[i])/(K1a+K1p);
      else
        for(int i=0; i<sim.sampleNr; i++) sim.c[ri].y[i]+=Vb*ca[i];
    }
 
  }
  /* 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); return(11);
    }
    int 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));
      return(12);
    }
    if(verbose>=0) {printf("%s saved.\n", simfile); fflush(stdout);}
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/