sim_3tcm.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 PET tissue time-radioactivity concentration curves (TACs)",
  "from arterial plasma (Ca) and blood (Cb) TACs, based on three-tissue",
  "compartmental model, where the compartments are in series (default):",
  " ",
  "  ____    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)                   ",
  "  Cvb(T) = Cab(T) - dCt(t)/dt / f                 ",
  "  Cpet(T)= Vb*fA*Cab(T) + Vb*(1-fA)*Cvb(T) + (1-Vb)*Ct(T) ",
  " ",
  ", or, optionally, three-tissue compartmental model, where the 2nd and 3rd",
  "tissue compartments are parallel, often used to represent specific and",
  "non-specific binding:",
  " ",
  "  ____    K1   ____    k3   ____   ",
  " | Ca | ----> | C1 | ----> | C2 |  ",
  " |____| <---- |____| <---- |____|  ",
  "          k2           k4          ",
  "                | ^                ",
  "             k5 | | k6             ",
  "                v |                ",
  "               ____                ",
  "              | C3 |               ",
  "              |____|               ",
  " ",
  "  dC1(t)/dt = K1*Ca(T) - (k2+k3+k5)*C1(T) + k4*C2(T) + k6*C3(T) ",
  "  dC2(t)/dt = k3*C1(T) - k4*C2(T)  ",
  "  dC3(t)/dt = k5*C1(T) - k6*C3(T)  ",
  "  Ct(T) = C1(T) + C2(T) + C3(T)    ",
  "  Cvb(T) = Cab(T) - dCt(t)/dt / f                 ",
  "  Cpet(T)= Vb*fA*Cab(T) + Vb*(1-fA)*Cvb(T) + (1-Vb)*Ct(T) ",
  " ",
  "Usage: @P [options] parfile [plasmafile bloodfile simfile]",
  " ",
  "Options:",
  " -paral[lel]",
  "     Model with parallel compartments C2 and C3 is applied, overriding",
  "     possible model settings in parameter file.",
  " -ser[ies]",
  "     Model with compartments C1, C2, and C3 in series is applied, overriding",
  "     possible model settings in parameter file.",
  " -vvm=<1|2>",
  "     Vascular volume is modelled as (1, default)",
  "     Cpet(T)= Vb*fA*Cab(T) + Vb*(1-fA)*Cvb(T) + (1-Vb)*Ct(T)",
  "     or as (2)",
  "     Cpet(T)= Vb*fA*Cab(T) + Vb*(1-fA)*Cvb(T) + Ct(T).",
  " -sub | -nosub",
  "     TACs of sub-compartments (C1, C2 and C3) are written (-sub)",
  "     or not written (-nosub, default) into the simfile.",
  " -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:",
  "K1, k2 (or K1/k2), k3, k4 (or k3/k4), k5, k6 (or k5/k6), and Vb.",
  "Additionally, values for perfusion (f), arterial fraction of Vb (fA), and",
  "delay time (dT) can be set; otherwise it is assumed that f>>K1, Cvb=Cab,",
  "fA=1, and dT=0.",
  "If parameter file does not contain units, then per min and per mL units",
  "are assumed, except sec for dT.",
  "For accurate results, plasma TAC should have very short sampling intervals.",
  "To reduce the model, k5 or k3 can be set to 0.",
  "When blood file is not used (Vb=0), 'none' can be entered in its place.",
  " ",
  "See also: simdisp, p2t_di, sim_rtcm, tacadd, simframe, tacformat, dft2img",
  " ",
  "Keywords: TAC, 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 main(int argc, char **argv)
{
  int    ai, help=0, version=0, verbose=1;
  char   plafile[FILENAME_MAX], blofile[FILENAME_MAX],
         simfile[FILENAME_MAX], parfile[FILENAME_MAX];
  char  *cptr;
  unsigned int model=0; // unknown
  unsigned int parNr=0;
  int    save_only_total=1;
  int    parallel=0;
  int    vvm=0; // Vascular volume model
  PAR    par;
  int    ret;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  plafile[0]=blofile[0]=simfile[0]=parfile[0]=(char)0;
  parInit(&par);
  //model=modelCodeIndex("SER3TCM");
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(!*cptr) continue;
    if(strncasecmp(cptr, "NOSUB", 5)==0) {
      save_only_total=1; continue;
    } else if(strncasecmp(cptr, "SUB", 3)==0) {
      save_only_total=0; continue;
    } else if(strncasecmp(cptr, "PARALLEL", 5)==0) {
      parallel=1; model=modelCodeIndex("PAR3TCM"); continue;
    } else if(strncasecmp(cptr, "SERIES", 3)==0) {
      parallel=0; model=modelCodeIndex("SER3TCM"); continue;
    } else if(strncasecmp(cptr, "VVM=", 4)==0) {
      if(atoiCheck(cptr+4, &vvm)==0 && vvm>=1 && vvm<=2) {vvm--; 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(plafile, 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(plafile[0] && !simfile[0]) {
    fprintf(stderr, "Error: missing filename.\n");
    return(1);
  }
  if(strcasecmp(blofile, "NONE")==0) strcpy(blofile, plafile);
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    if(model!=0) printf("model := %s\n", modelCode(model));
    printf("parfile := %s\n", parfile);
    printf("plafile := %s\n", plafile);
    printf("blofile := %s\n", blofile);
    printf("simfile := %s\n", simfile);
    printf("save_only_total := %d\n", save_only_total);
    printf("model := %s\n", modelCode(model));
    if(model!=0) printf("parallel := %d\n", parallel);
    printf("vvm := %d\n", 1+vvm);
  }
  
  /*
   *  Make template parameter file, if no other filenames were given,
   *  and then quit
   */
  if(!plafile[0]) {
    if(model==0) {parallel=0; model=modelCodeIndex("SER3TCM");}
    parNr=10;
    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;
    if(parallel) {
      strcpy(par.n[0].name, "K1"); par.n[0].unit=UNIT_ML_PER_ML_MIN;
      strcpy(par.n[1].name, "K1/k2"); par.n[1].unit=UNIT_ML_PER_ML;
      strcpy(par.n[2].name, "k3"); par.n[2].unit=UNIT_PER_MIN;
      strcpy(par.n[3].name, "k3/k4"); par.n[3].unit=UNIT_UNITLESS;
      strcpy(par.n[4].name, "k5"); par.n[4].unit=UNIT_PER_MIN;
      strcpy(par.n[5].name, "k5/k6"); par.n[5].unit=UNIT_UNITLESS;
    } else {
      strcpy(par.n[0].name, "K1"); par.n[0].unit=UNIT_ML_PER_ML_MIN;
      strcpy(par.n[1].name, "k2"); par.n[1].unit=UNIT_PER_MIN;
      strcpy(par.n[2].name, "k3"); par.n[2].unit=UNIT_PER_MIN;
      strcpy(par.n[3].name, "k4"); par.n[3].unit=UNIT_PER_MIN;
      strcpy(par.n[4].name, "k5"); par.n[4].unit=UNIT_PER_MIN;
      strcpy(par.n[5].name, "k6"); par.n[5].unit=UNIT_PER_MIN;
    }
    strcpy(par.n[6].name, "Vb"); par.n[6].unit=UNIT_PERCENTAGE;
    strcpy(par.n[7].name, "f"); par.n[7].unit=UNIT_ML_PER_ML_MIN;
    strcpy(par.n[8].name, "fA"); par.n[8].unit=UNIT_PERCENTAGE;
    strcpy(par.n[9].name, "dT"); par.n[9].unit=UNIT_SEC;
    for(int i=0; i<par.tacNr; i++) {
      sprintf(par.r[i].name, "tac%d", 1+i);
      par.r[i].model=model;
      //par.r[i].start=0.0;
      //par.r[i].end=90.0;
      par.r[i].p[0]=0.3;
      par.r[i].p[1]=0.8;
      par.r[i].p[6]=4.0;
      par.r[i].p[7]=0.8;
      par.r[i].p[8]=30.;
      par.r[i].p[9]=0.0;
    }
    if(parallel) {
      int i=0;
      par.r[i].p[2]=0.0;
      par.r[i].p[3]=0.0;
      par.r[i].p[4]=0.5;
      par.r[i].p[5]=1.0;
      i=1;
      par.r[i].p[2]=0.15;
      par.r[i].p[3]=1.5;
      par.r[i].p[4]=0.5;
      par.r[i].p[5]=1.0;
      i=2;
      par.r[i].p[2]=0.15;
      par.r[i].p[3]=3.0;
      par.r[i].p[4]=0.5;
      par.r[i].p[5]=1.0;
    } else {
      int i=0;
      par.r[i].p[2]=0.20;
      par.r[i].p[3]=0.15;
      par.r[i].p[4]=0.01;
      par.r[i].p[5]=0.0;
      i=1;
      par.r[i].p[2]=0.20;
      par.r[i].p[3]=0.15;
      par.r[i].p[4]=0.05;
      par.r[i].p[5]=0.0;
      i=2;
      par.r[i].p[2]=0.20;
      par.r[i].p[3]=0.15;
      par.r[i].p[4]=0.10;
      par.r[i].p[5]=0.01;
    }
    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); 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);
  if(model!=0) {
    /* If user gave model, then use it */
    for(int i=0; i<par.tacNr; i++) par.r[i].model=model;
  } else {
    /* otherwise, check the model in parameter file */
    for(int i=0; i<par.tacNr; i++) {
      /* If parameter file does not contain model, then use default */
      if(par.r[i].model==0) { 
        par.r[i].model=modelCodeIndex("SER3TCM");
      } else {
        /* If parameter file contains model, then check that it is ok here */
        if(par.r[i].model==modelCodeIndex("SER3TCM")) continue;
        if(par.r[i].model==modelCodeIndex("PAR3TCM")) continue;
        if(par.r[i].model==modelCodeIndex("2TCM")) continue;
        if(par.r[i].model==modelCodeIndex("1TCM")) continue;
        fprintf(stderr, "Error: model %s not supported.\n", modelCode(par.r[i].model));
        parFree(&par); fflush(stderr); return(3);
      }
    }
  }
  /* Check the validity of model and its parameters */
  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);
    /* check that we got at least as many parameters as the model has, minus
       the three optional parameters */
    if(modelParNr(par.r[i].model)>4+(unsigned int)par.parNr) {
      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, "K1")<0) ret++;
    /*
    if(parFindParameter(&par, "k2")<0 && parFindParameter(&par, "K1/k2")<0) ret++;
    if(parFindParameter(&par, "k3")<0) ret++;
    if(parFindParameter(&par, "k4")<0 && parFindParameter(&par, "k3/k4")<0) ret++;
    if(parFindParameter(&par, "k5")<0) ret++;
    if(parFindParameter(&par, "k6")<0 && parFindParameter(&par, "k5/k6")<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 TACs\n");
  TAC input; tacInit(&input);
  ret=tacReadModelingInput(plafile, blofile, 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"); fflush(stderr);
  }
 
 
  /*
   *  Check parameter and data units
   */
  {
    if(verbose>1) {fprintf(stdout, "checking parameters and data units\n"); fflush(stdout);}
    int i, punit;
    /* Vb */
    if((i=parFindParameter(&par, "Vb"))>=0) {
      if(verbose>2) {printf("  Vb\n"); fflush(stdout);}
      punit=par.n[i].unit;
      if(punit==UNIT_UNKNOWN) {
        double pmax=0.0; (void)parValueRange(&par, i, NULL, &pmax);
        if(pmax>1.0) punit=UNIT_PERCENTAGE; else punit=UNIT_UNITLESS;
      }
      if(punit==UNIT_PERCENTAGE) {
        if(verbose>1) printf("Note: converting Vb from percentage to fractions.\n");
        par.n[i].unit=UNIT_UNITLESS;
        for(int j=0; j<par.tacNr; j++) par.r[j].p[i]*=0.01;
      }
    }
    /* fA */
    if((i=parFindParameter(&par, "fA"))>=0) {
      if(verbose>2) {printf("  fA\n"); fflush(stdout);}
      punit=par.n[i].unit;
      if(punit==UNIT_UNKNOWN) {
        double pmax=0.0; (void)parValueRange(&par, i, NULL, &pmax);
        if(pmax>1.0) punit=UNIT_PERCENTAGE; else punit=UNIT_UNITLESS;
      }
      if(punit==UNIT_PERCENTAGE) {
        if(verbose>1) printf("Note: converting fA from percentage to fractions.\n");
        par.n[i].unit=UNIT_UNITLESS;
        for(int j=0; j<par.tacNr; j++) par.r[j].p[i]*=0.01;
      }
    }
    /* Perfusion */
    i=parFindParameter(&par, "f"); if(i<0) i=parFindParameter(&par, "Flow");
    if(i>=0) {
      if(verbose>2) {printf("  f\n"); fflush(stdout);}
      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.;
      }
    }
    /* dT */
    if((i=parFindParameter(&par, "dT"))>=0) {
      if(verbose>2) {printf("  dT\n"); fflush(stdout);}
      punit=par.n[i].unit;
      if(verbose>4) printf("  units: dT=%s x=%s\n", unitName(punit), unitName(input.tunit));
      if(punit==UNIT_UNKNOWN) punit=UNIT_SEC;
      if(input.tunit==UNIT_MIN && punit==UNIT_SEC) {
        if(verbose>1) printf("Note: converting dT from sec to min.\n");
        par.n[i].unit=UNIT_MIN;
        for(int j=0; j<par.tacNr; j++) par.r[j].p[i]/=60.;
      } else if(input.tunit==UNIT_SEC && punit==UNIT_MIN) {
        if(verbose>1) printf("Note: converting dT from min to sec.\n");
        par.n[i].unit=UNIT_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) {
        if(verbose>2) {printf("  %s\n", rcname); fflush(stdout);}
        punit=par.n[i].unit;
        if(verbose>5) 
          for(int j=0; j<par.tacNr; j++)
            printf("  %s %s=%g %s\n", par.r[j].name, rcname, par.r[j].p[i], unitName(punit));
        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.;
        }
      }
    }
    if(verbose>3) parWrite(&par, stdout, PAR_FORMAT_TSV_UK, 1, NULL);
  }
 
 
  /*
   *  Allocate space for simulated data
   */
  if(verbose>1) fprintf(stdout, "allocating space for simulated TACs\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];
 
 
  /*
   *  Simulation
   */
  if(verbose>1) printf("simulating\n");
  double *px=input.x;
  double *ppy=input.c[0].y;
  double *pby=input.c[0].y; if(input.tacNr>1) pby=input.c[1].y;
  for(int i=0; i<par.tacNr; i++) {
    if(verbose>2) printf("simulating %s\n", sim.c[i].name);
    /* Check whether this tac model is parallel */
    if(par.r[i].model==modelCodeIndex("PAR3TCM")) parallel=1; else parallel=0;
    /* Create pointers for simulated data, but set those only after possible reallocation */
    double *ppet=sim.c[i].y; // set again if reallocated
    double *pc1=NULL;
    double *pc2=NULL;
    double *pc3=NULL;
    double *pcab=NULL;
    double *pcvb=NULL;
    /* Get parameters */
    double K1, k2, k3, k4, k5, k6, Vb, Flow, fA, dT;
    K1=k2=k3=k4=k5=k6=Vb=Flow=fA=dT=nan("");
    K1=parGetParameter(&par, "K1", i);
    k2=parGetParameter(&par, "k2", i);
    if(isnan(k2)) {
      double f=parGetParameter(&par, "K1/k2", i);
      if(isnormal(f)) k2=K1/f; // not if zero, subnormal, infinite, or NaN
    }
    k3=parGetParameter(&par, "k3", i);
    k4=parGetParameter(&par, "k4", i);
    if(isnan(k4)) {
      double f=parGetParameter(&par, "k3/k4", i);
      if(isnormal(f)) k4=k3/f; // not if zero, subnormal, infinite, or NaN
    }
    k5=parGetParameter(&par, "k5", i);
    k6=parGetParameter(&par, "k6", i);
    if(isnan(k6)) {
      double f=parGetParameter(&par, "k5/k6", i);
      if(isnormal(f)) k6=k5/f; // not if zero, subnormal, infinite, or NaN
    }
    Vb=parGetParameter(&par, "Vb", i);
    Flow=parGetParameter(&par, "f", i);
    if(isnan(Flow)) Flow=parGetParameter(&par, "Flow", i);
    fA=parGetParameter(&par, "fA", i);
    dT=parGetParameter(&par, "dT", i);
    /* Verify and fix parameters */
    if(!(K1>=0.0) || k2<0 || k3<0 || k4<0 || k5<0 || k6<0) {
      fprintf(stderr, "Error: invalid rate constant.\n");
      tacFree(&sim); tacFree(&input); parFree(&par); return(5);
    }
    if(K1==0.0) {k2=k3=k4=k5=k6=0.0;}
    if(k3==0.0) {k4=0.0; if(!parallel) k5=k6=0.0;}
    if(k5==0.0) k6=0.0;
    if(isnan(Vb) || Vb==0.0) {Flow=fA=nan("");}
    if(isnan(Flow) || Flow==0.0 || isnan(fA) || fA==0.0 || fA==1.0) {Flow=fA=nan("");}
    if(Vb<0.0 || Vb>1.0) {
      fprintf(stderr, "Error: invalid vascular volume.\n");
      tacFree(&sim); tacFree(&input); parFree(&par); return(5);
    }
    if(Flow<K1) {
      fprintf(stderr, "Error: invalid perfusion value (f<K1).\n");
      tacFree(&sim); tacFree(&input); parFree(&par); return(5);
    }
    if(fA<0.0 || fA>1.0) {
      fprintf(stderr, "Error: invalid fA value.\n");
      tacFree(&sim); tacFree(&input); parFree(&par); return(5);
    }
    if(verbose>3) {
      printf("K1 := %g\n", K1);
      printf("k2 := %g\n", k2);
      printf("k3 := %g\n", k3);
      printf("k4 := %g\n", k4);
      printf("k5 := %g\n", k5);
      printf("k6 := %g\n", k6);
      printf("Vb := %g\n", Vb);
      printf("f := %g\n", Flow);
      printf("fA := %g\n", fA);
      printf("dT := %g\n", dT);
    }
    /* Add place for compartment TACs in simulated data, as needed */
    int currentTacNr=1;
    if(save_only_total==0) {
      int n=0; 
      if(Vb>0.0 || k3>0.0 || k5>0.0) n++; 
      if(k3>0.0) n++;
      if(k5>0.0) n++;
      if(Vb>0.0) n++;
      if(Vb>0.0 && fA<1.0 && Flow>0.0) n++;
      if(n>0 && tacAllocateMore(&sim, n)!=0) {
        fprintf(stderr, "Error: cannot allocate memory for simulations.\n");
        tacFree(&input); tacFree(&sim); parFree(&par); return(4);
      }
      ppet=sim.c[i].y;
      if(n>0) {
        strcpy(sim.c[sim.tacNr].name, sim.c[i].name);
        strcat(sim.c[sim.tacNr].name, "_C1");
        pc1=sim.c[sim.tacNr].y;
        sim.tacNr++; n--; currentTacNr++;
      } 
      if(n>0 && k3>0.0) {
        strcpy(sim.c[sim.tacNr].name, sim.c[i].name);
        strcat(sim.c[sim.tacNr].name, "_C2");
        pc2=sim.c[sim.tacNr].y;
        sim.tacNr++; n--; currentTacNr++;
      }
      if(n>0 && k5>0.0) {
        strcpy(sim.c[sim.tacNr].name, sim.c[i].name);
        strcat(sim.c[sim.tacNr].name, "_C3");
        pc3=sim.c[sim.tacNr].y;
        sim.tacNr++; n--; currentTacNr++;
      }
      /* Vascular TACs */
      if(n>0 && Vb>0.0) {
        strcpy(sim.c[sim.tacNr].name, sim.c[i].name);
        strcat(sim.c[sim.tacNr].name, "_Cab");
        pcab=sim.c[sim.tacNr].y;
        sim.tacNr++; n--; currentTacNr++;
      } 
      if(n>0 && Vb>0.0 && fA<1.0 && Flow>0.0) {
        strcpy(sim.c[sim.tacNr].name, sim.c[i].name);
        strcat(sim.c[sim.tacNr].name, "_Cvb");
        pcvb=sim.c[sim.tacNr].y;
        sim.tacNr++; n--; currentTacNr++;
      } 
    }
    /* Simulate */
    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;
    if(Vb==1.0) {
      for(int i=0; i<sim.sampleNr; i++) ppet[i]=pby[i];
      if(isnan(fA)) fA=1.0;
      if(pcab!=NULL) for(int i=0; i<sim.sampleNr; i++) pcab[i]=fA*pby[i];
      if(fA<1.0 && pcab!=NULL) for(int i=0; i<sim.sampleNr; i++) pcvb[i]=(1.0-fA)*pby[i];
      ret=0;
    } else if(Vb>0.0) {
      if(isnan(Flow)) Flow=0.0;
      if(isnan(fA)) fA=1.0;
      if(parallel==0)
        ret=simC3vs(px, ppy, pby, sim.sampleNr, K1, k2, k3, k4, k5, k6, Flow, Vb, fA, vvm,
                    ppet, pc1, pc2, pc3, pcab, pcvb);
      else
        ret=simC3vp(px, ppy, pby, sim.sampleNr, K1, k2, k3, k4, k5, k6, Flow, Vb, fA, vvm,
                    ppet, pc1, pc2, pc3, pcab, pcvb);
    } else {
      if(parallel==0)
        ret=simC3s(px, ppy, sim.sampleNr, K1, k2, k3, k4, k5, k6, 
                   ppet, pc1, pc2, pc3);
      else
        ret=simC3p(px, ppy, sim.sampleNr, K1, k2, k3, k4, k5, k6,
                   ppet, pc1, pc2, pc3);
    }
    if(ret) {
      fprintf(stderr, "Error: invalid data for simulation.\n");
      if(verbose>1) printf("sim_return_code := %d\n", ret);
      tacFree(&sim); tacFree(&input); parFree(&par); return(6);
    }
    /* Simulate delay time */
    if(!isnan(dT) && fabs(dT)>1.0E-20) {
      for(int j=0; j<currentTacNr; j++) {
        if(verbose>4) {
          printf("  delay for %s[%d/%d]\n", sim.c[i].name, 1+j, currentTacNr); fflush(stdout);}
        if((ret=tacDelay(&sim, dT, i+j, &status)) != TPCERROR_OK) break;
      }
      if(ret) {
        fprintf(stderr, "Error: invalid delay time for simulation.\n");
        if(verbose>2) printf("sim_return_code := %d\n", ret);
        tacFree(&sim); tacFree(&input); parFree(&par); return(7);
      }
    }
  } // next TAC
 
  /* 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);
  }
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/