p2t_di.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Simulation of PET tissue time-radioactivity concentration curve (TTAC) using",
  "dual-input four-tissue compartmental models. Input consists of parent tracer",
  "and metabolite concentration curves in arterial plasma (Cap and Cam) and",
  "total radioactivity concentration in arterial blood (Cb).",
  " ",
  "Tissue compartments for parent tracer are in series [1] (by default):",
  " ",
  "  _____    K1   ____    k3   ____    k5   ____    k7  ",
  " | Cap | ----> | C1 | ----> | C2 | ----> | C3 | ----> ",
  " |_____| <---- |____| <---- |____| <---- |____|       ",
  "           k2     |     k4           k6               ",
  "                km|                                   ",
  "                  v                                   ",
  "  _____   K1m   ____                                  ",
  " | Cam | ----> | C4 |                                 ",
  " |_____| <---- |____|                                 ",
  "          k2m                                         ",
  " ",
  "  dC1(t)/dt = K1*Cap(T) - (k2+k3+km)*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+k7)*C3(T)                ",
  "  dC4(t)/dt = K1m*Cam(T) + km*C1(T) - k2m*C4(T)       ",
  "  Ct(T) = C1(T) + C2(T) + C3(T) + C4(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, the 2nd and 3rd tissue compartments are parallel [2],",
  "often used to represent specific and non-specific binding:",
  " ",
  "                ____               ",
  "               | C3 |   k7         ",
  "               |____| ---->        ",
  "                 ^ |               ",
  "              k5 | | k6            ",
  "                 | v               ",
  "  _____    K1   ____    k3   ____  ",
  " | Cap | ----> | C1 | ----> | C2 | ",
  " |_____| <---- |____| <---- |____| ",
  "           k2     |     k4         ",
  "                km|                ",
  "                  v                ",
  "  _____   K1m   ____               ",
  " | Cam | ----> | C4 |              ",
  " |_____| <---- |____|              ",
  "          k2m                      ",
  " ",
  "  dC1(t)/dt = K1*Ca(T) - (k2+k3+k5+km)*C1(T) + k4*C2(T) + k6*C3(T) ",
  "  dC2(t)/dt = k3*C1(T) - k4*C2(T)                ",
  "  dC3(t)/dt = k5*C2(T) - (k6+k7)*C3(T)           ",
  "  dC4(t)/dt = K1m*Cam(T) + km*C1(T) - k2m*C4(T)  ",
  "  Ct(T) = C1(T) + C2(T) + C3(T) + C4(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] parentfile metabolitefile bloodfile"
  " K1 k2 k3 k4 k5 k6 k7 km K1m k2m Vb simfile",
  " ",
  "Options:",
  " -paral[lel]",
  "     Model with parallel compartments C2 and C3 is applied.",
  " -ser[ies]",
  "     Model with compartments C1, C2, and C3 in series is applied (default).",
  " -sub | -nosub",
  "     TACs of sub-compartments (C1, C2 and C3) are written (-sub)",
  "     or not written (-nosub, default) into the output file.",
  " -f=<Perfusion (ml/(min*ml) or ml/(sec*ml))>",
  "     Difference between concentrations in venous and arterial blood can",
  "     be simulated if tissue perfusion (f>K1) is specified with this option",
  "     and arterial fraction of vascular volume is set with option -fA;",
  "     by default it is assumed that venous and arterial activities are",
  "     the same (f>>K1).",
  " -fA=<Arterial fraction of vascular volume (%)>",
  "     Difference between concentrations in venous and arterial blood can",
  "     be simulated if arterial fraction of vascular volume is specified with",
  "     this option and tissue perfusion is set with option -f;",
  "     by default it is assumed that Vb consists of only arterial blood.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "If the times in plasma file are in seconds, the units of rate constants",
  "(k's) and blood flow (f) must also be specified as 1/sec.",
  "For accurate results, plasma TAC should have very short sampling intervals.",
  "To reduce the model, k7, k5, and k3 can be set to 0.",
  " ",
  "Simulated TACs are written in ASCII format with columns:",
  "  1) Sample time",
  "  2) Total tissue activity concentration (Cpet)",
  "  3) Activity concentration in 1st tissue compartment, (1-Vb)*C1 (optional)",
  "  4) Activity concentration in 2nd tissue compartment, (1-Vb)*C2 (optional)",
  "  5) Activity concentration in 3rd tissue compartment, (1-Vb)*C3 (optional)",
  "  6) Activity concentration in 4th tissue compartment, (1-Vb)*C4 (optional)",
  "  7) Arterial contribution to tissue activity, Vb*fA*Cab (optional)",
  "  8) Venous contribution to tissue activity, Vb*(1-fA)*Cvb (optional)",
  " ",
  "References:",
  "1. TPCMOD0001 Appendix B.",
  "2. TPCMOD0001 Appendix C.",
  " ",
  "See also: sim_3tcm, tacadd, tacren, simframe, tacformat, dft2img",
  " ",
  "Keywords: TAC, simulation, modelling, compartmental model, dual-input",
  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;
  int        fi, ri, ret, save_only_total=1;
  int        parallel=0;
  DFT        plasma, sim;
  double     k1, k2, k3, k4, k5, k6, k7, km, k1m, k2m, f, Vb, fA;
  double    *t, *cap, *cam, *cab, *cpet, *ct1, *ct2, *ct3, *ct4, *ctab, *ctvb;
  char       pfile[FILENAME_MAX], mfile[FILENAME_MAX], bfile[FILENAME_MAX],
             sfile[FILENAME_MAX], *cptr, tmp[FILENAME_MAX+128];
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  pfile[0]=mfile[0]=bfile[0]=sfile[0]=(char)0;
  k1=k2=k3=k4=k5=k6=k7=km=k1m=k2m=Vb=-1.0;
  f=0.0; fA=1.0;
  dftInit(&plasma); dftInit(&sim);
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    cptr=argv[ai]+1;
    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; continue;
    } else if(strncasecmp(cptr, "SERIES", 3)==0) {
      parallel=0; continue;
    } else if(strncasecmp(cptr, "F=", 2)==0) {
      ret=atof_with_check(cptr+2, &f); if(ret==0 && f>=0.0) continue;
    } else if(strncasecmp(cptr, "FA=", 3)==0) {
      ret=atof_with_check(cptr+3, &fA);
      if(ret==0 && fA>=0.0 && fA<100.0) {fA*=0.01; continue;}
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break;
 
  /* 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 */
  for(; ai<argc; ai++) {
    if(!pfile[0]) {
      strcpy(pfile, argv[ai]); continue;
    } else if(!mfile[0]) {
      strcpy(mfile, argv[ai]); continue;
    } else if(!bfile[0]) {
      strcpy(bfile, argv[ai]); continue;
    } else if(k1<0.0) {
      k1=atof_dpi(argv[ai]); if(k1<0.0) k1=0.0; continue;
    } else if(k2<0.0) {
      k2=atof_dpi(argv[ai]); if(k2<0.0) k2=0.0; continue;
    } else if(k3<0.0) {
      k3=atof_dpi(argv[ai]); if(k3<0.0) k3=0.0; continue;
    } else if(k4<0.0) {
      k4=atof_dpi(argv[ai]); if(k4<0.0) k4=0.0; continue;
    } else if(k5<0.0) {
      k5=atof_dpi(argv[ai]); if(k5<0.0) k5=0.0; continue;
    } else if(k6<0.0) {
      k6=atof_dpi(argv[ai]); if(k6<0.0) k6=0.0; continue;
    } else if(k7<0.0) {
      k7=atof_dpi(argv[ai]); if(k7<0.0) k7=0.0; continue;
    } else if(km<0.0) {
      km=atof_dpi(argv[ai]); if(km<0.0) km=0.0; continue;
    } else if(k1m<0.0) {
      k1m=atof_dpi(argv[ai]); if(k1m<0.0) k1m=0.0; continue;
    } else if(k2m<0.0) {
      k2m=atof_dpi(argv[ai]); if(k2m<0.0) k2m=0.0; continue;
    } else if(Vb<0.0) {
      ret=atof_with_check(argv[ai], &Vb);
      if(ret==0 && Vb>=0.0 && Vb<100.0) {Vb*=0.01; continue;}
    } else if(!sfile[0]) {
      strcpy(sfile, argv[ai]); continue;
    } else {
      fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
      return(1);
    }
  }
 
  /* Is something missing? */
  if(!sfile[0] || Vb<0.0) {tpcPrintUsage(argv[0], info, stdout); return(1);}
  /* Or otherwise wrong? */
  if(f>0.0 && k1>0.0 && f<k1) {
    fprintf(stderr, "Error: f cannot be lower than K1!\n"); return(1);}
  if(f>0.0 && k1m>0.0 && f<k1m) {
    fprintf(stderr, "Error: f cannot be lower than K1m!\n"); return(1);}
  if(Vb<=0.0 && (fA<1.0 || f>0.0))
    fprintf(stderr, "Warning: options -fA and -f are ignored when Vb=0%%\n");
  else if(f>0.0 && fA==1.0)
    fprintf(stderr, "Warning: option -f is ignored when fA=100%%\n");
  if(k1>1.0E+08 || k2>1.0E+08 || k3>1.0E+08 || k4>1.0E+08 || k5>1.0E+08 || 
     k6>1.0E+08 || k7>1.0E+08 || km>1.0E+08 || k1m>1.0E+08 || k2m>1.0E+08) {
    fprintf(stderr, "Error: too high rate constant.\n"); return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("pfile := %s\n", pfile);
    printf("mfile := %s\n", mfile);
    printf("bfile := %s\n", bfile);
    printf("sfile := %s\n", sfile);
    printf("rate_constants := %g %g %g %g %g %g %g %g %g %g\n",
      k1, k2, k3, k4, k5, k6, k7, km, k1m, k2m);
    printf("Vb := %g\n", Vb);
    printf("fA := %g\n", fA);
    printf("f := %g\n", f);
    printf("save_only_total := %d\n", save_only_total);
    printf("parallel := %d\n", parallel);
  }
 
 
  /*
   *  Read plasma parent data
   */
  if(verbose>1) printf("reading plasma parent TAC\n");
  if(dftRead(pfile, &plasma)) {
    fprintf(stderr, "Error in reading '%s': %s\n", pfile, dfterrmsg);
    dftEmpty(&plasma); return(2);
  }
  if(plasma.frameNr<3) {
    fprintf(stderr, "Error: too few samples in plasma data.\n");
    dftEmpty(&plasma); return(2);
  }
  if(plasma.voiNr>1) {
    fprintf(stderr, "Error: plasma data contains more than one curve.\n");
    dftEmpty(&plasma); return(2);
  }
  if(plasma.timeunit==TUNIT_UNKNOWN) {
    plasma.timeunit=TUNIT_MIN;
    if(verbose>=0) printf("assuming that time unit is min\n");
  }
  strcpy(plasma.voi[0].voiname, "Parent");
  strcpy(plasma.voi[0].name, plasma.voi[0].voiname);
 
 
  /*
   *  Read plasma metabolite data; interpolate it to the parent times
   */
  if(verbose>1) printf("reading plasma metabolite TAC\n");
  if(dftRead(mfile, &sim)) {
    fprintf(stderr, "Error in reading '%s': %s\n", mfile, dfterrmsg);
    dftEmpty(&plasma); dftEmpty(&sim); return(3);
  }
  if(sim.frameNr<3) {
    fprintf(stderr, "Error: too few samples in metabolite data.\n");
    dftEmpty(&plasma); dftEmpty(&sim); return(3);
  }
  if(sim.voiNr>1) {
    fprintf(stderr, "Error: metabolite data contains more than one curve.\n");
    dftEmpty(&plasma); dftEmpty(&sim); return(3);
  }
  /* Check that blood has been measured for long enough */
  if(sim.x[sim.frameNr-1] < 0.85*plasma.x[plasma.frameNr-1]) {
    fprintf(stderr, "Warning: metabolite TAC is shorter than parent TAC.\n");
  }
  /* Allocate memory for interpolated metabolite data */
  if(dftAddmem(&plasma, 1)) {
    fprintf(stderr, "Error: cannot allocate memory for metabolite TAC.\n");
    dftEmpty(&plasma); dftEmpty(&sim); return(3);
  }
  strcpy(plasma.voi[1].voiname, "Metab");
  strcpy(plasma.voi[1].name, plasma.voi[1].voiname);
  /* Interpolate */
  ret=interpolate(sim.x, sim.voi[0].y, sim.frameNr, plasma.x, plasma.voi[1].y,
      plasma.voi[1].y2, NULL, plasma.frameNr);
  dftEmpty(&sim);
  if(ret) {
    fprintf(stderr, "Error: cannot interpolate metabolite data.\n");
    dftEmpty(&plasma); return(3);
  }
  plasma.voiNr++;
  if(verbose>20) dftPrint(&plasma);
  /* Verify that metabolite TAC is > 0 */
  if(plasma.voi[1].y2[plasma.frameNr-1]<=0.0) {
    fprintf(stderr, "Warning: plasma metabolite concentration <= 0.\n");
  }
 
 
  /*
   *  Read blood data; interpolate it to the plasma times
   */
  if(verbose>1) printf("reading blood\n");
  if(dftRead(bfile, &sim)) {
    fprintf(stderr, "Error in reading '%s': %s\n", bfile, dfterrmsg);
    dftEmpty(&plasma); dftEmpty(&sim); return(3);
  }
  if(sim.frameNr<3) {
    fprintf(stderr, "Error: too few samples in blood data.\n");
    dftEmpty(&plasma); dftEmpty(&sim); return(3);
  }
  if(sim.voiNr>1) {
    fprintf(stderr, "Error: blood data contains more than one curve.\n");
    dftEmpty(&plasma); dftEmpty(&sim); return(3);
  }
  /* Check that blood has been measured for long enough */
  if(sim.x[sim.frameNr-1] < 0.67*plasma.x[plasma.frameNr-1]) {
    fprintf(stderr, "Warning: blood TAC is much shorter than plasma TAC.\n");
  }
  /* Allocate memory for interpolated blood data */
  if(dftAddmem(&plasma, 1)) {
    fprintf(stderr, "Error: cannot allocate memory for blood TAC.\n");
    dftEmpty(&plasma); dftEmpty(&sim); return(3);
  }
  strcpy(plasma.voi[2].voiname, "Blood");
  strcpy(plasma.voi[2].name, plasma.voi[2].voiname);
  /* Interpolate */
  ret=interpolate(sim.x, sim.voi[0].y, sim.frameNr, plasma.x, plasma.voi[2].y,
      plasma.voi[2].y2, NULL, plasma.frameNr);
  dftEmpty(&sim);
  if(ret) {
    fprintf(stderr, "Error: cannot interpolate blood data.\n");
    dftEmpty(&plasma); return(3);
  }
  plasma.voiNr++;
  if(verbose>8) dftPrint(&plasma);
  /* Verify that blood TAC is > 0 */
  if(plasma.voi[2].y2[plasma.frameNr-1]<=0.0) {
    fprintf(stderr, "Warning: blood concentration <= 0.\n");
  }
 
 
  if(verbose>5) {
    strcpy(tmp, "simulation_input.dat");
    printf("Saving input in %s for testing.\n", tmp);
    if(dftWrite(&plasma, tmp)) {
      fprintf(stderr, "Error in writing '%s': %s\n", tmp, dfterrmsg);
    }
  }
  
 
  /*
   *  Allocate memory for simulated tissue data.
   */
  if(verbose>1) printf("allocating memory\n");
  /* allocate */
  ret=dftSetmem(&sim, plasma.frameNr, 7);
  if(ret) {
    fprintf(stderr, 
            "Error (%d): cannot allocate memory for simulated TACs.\n", ret);
    dftEmpty(&plasma); return(4);
  }
  sim.frameNr=plasma.frameNr; sim.voiNr=7;
  /* Copy times & header info */
  dftCopymainhdr(&plasma, &sim);
  for(fi=0; fi<sim.frameNr; fi++) {
    sim.x[fi]=plasma.x[fi];
    sim.x1[fi]=plasma.x1[fi]; sim.x2[fi]=plasma.x2[fi];
  }
  /* column headers */
  for(ri=0; ri<sim.voiNr; ri++) {
    sim.voi[ri].sw=0;
    switch(ri) {
      case 0: strcpy(sim.voi[ri].voiname, "Cpet");
        sim.voi[ri].size=100.0;
        break;
      case 1: strcpy(sim.voi[ri].voiname, "C1"); 
        sim.voi[ri].size=100.*(1.-Vb);
        break;
      case 2: strcpy(sim.voi[ri].voiname, "C2"); 
        sim.voi[ri].size=100.*(1.-Vb);
        break;
      case 3: strcpy(sim.voi[ri].voiname, "C3");
        sim.voi[ri].size=100.*(1.-Vb);
        break;
      case 4: strcpy(sim.voi[ri].voiname, "C4");
        sim.voi[ri].size=100.*(1.-Vb);
        break;
      case 5: strcpy(sim.voi[ri].voiname, "Cab");
        sim.voi[ri].size=100.0*fA*Vb;
        break;
      case 6: strcpy(sim.voi[ri].voiname, "Cvb");
        sim.voi[ri].size=100.0*(1.0-fA)*Vb;
        break;
    }
    strcpy(sim.voi[ri].hemisphere, ""); strcpy(sim.voi[ri].place, "");
    strcpy(sim.voi[ri].name, sim.voi[ri].voiname);
  }
 
  /*
   *  Simulate TACs
   */
  if(verbose>1) printf("simulating\n");
  t=sim.x; cap=plasma.voi[0].y; cam=plasma.voi[1].y; cab=plasma.voi[2].y;
  cpet=sim.voi[0].y; ct1=sim.voi[1].y; ct2=sim.voi[2].y; ct3=sim.voi[3].y;
  ct4=sim.voi[4].y; ctab=sim.voi[5].y; ctvb=sim.voi[6].y;
  if(parallel==0)
    ret=simC4DIvs(t, cap, cam, cab, sim.frameNr,
                  k1, k2, k3, k4, k5, k6, k7, km, k1m, k2m, f, Vb, fA,
                  cpet, ct1, ct2, ct3, ct4, ctab, ctvb, verbose-1);
  else
    ret=simC4DIvp(t, cap, cam, cab, sim.frameNr,
                  k1, k2, k3, k4, k5, k6, k7, km, k1m, k2m, f, Vb, fA,
                  cpet, ct1, ct2, ct3, ct4, ctab, ctvb, verbose-1);
  dftEmpty(&plasma);
  if(ret!=0) {
    fprintf(stderr, "Error (%d) in simulation.\n", ret);
    dftEmpty(&sim); return(8);
  }
 
 
  /*
   *  Save simulated TACs
   */
  if(verbose>1) printf("saving PET curves\n");
  DFT_NR_OF_DECIMALS=6;
  if(save_only_total) sim.voiNr=1;
  else if(Vb<=0.0) sim.voiNr=5;
  else sim.voiNr=7;
  /* Set comments */
  dftSetComments(&sim);
  sprintf(tmp, "# parentfile := %s\n", pfile); strcat(sim.comments, tmp);
  sprintf(tmp, "# metabolitefile := %s\n", mfile); strcat(sim.comments, tmp);
  sprintf(tmp, "# bloodfile := %s\n", bfile); strcat(sim.comments, tmp);
  strcpy(tmp, "# model := ");
  if(parallel==0) strcat(tmp, "C4DIVS\n"); else strcat(tmp, "C4DIVP\n");
  strcat(sim.comments, tmp);
  sprintf(tmp, "# K1 := %g\n", k1); strcat(sim.comments, tmp);
  sprintf(tmp, "# k2 := %g\n", k2); strcat(sim.comments, tmp);
  sprintf(tmp, "# k3 := %g\n", k3); strcat(sim.comments, tmp);
  sprintf(tmp, "# k4 := %g\n", k4); strcat(sim.comments, tmp);
  sprintf(tmp, "# k5 := %g\n", k5); strcat(sim.comments, tmp);
  sprintf(tmp, "# k6 := %g\n", k6); strcat(sim.comments, tmp);
  sprintf(tmp, "# k7 := %g\n", k7); strcat(sim.comments, tmp);
  sprintf(tmp, "# km := %g\n", km); strcat(sim.comments, tmp);
  sprintf(tmp, "# K1m := %g\n", k1m); strcat(sim.comments, tmp);
  sprintf(tmp, "# k2m := %g\n", k2m); strcat(sim.comments, tmp);
  if(Vb>0.0) {
    sprintf(tmp, "# f := %g\n", f); strcat(sim.comments, tmp);
    sprintf(tmp, "# Vb := %g [%%]\n", 100.0*Vb); strcat(sim.comments, tmp);
    sprintf(tmp, "# fA := %g [%%]\n", 100.0*fA); strcat(sim.comments, tmp);
  }
  /* Set file format to PMOD, if extension is .tac */
  if(!strcasecmp(filenameGetExtension(sfile), ".tac"))
    sim._type=DFT_FORMAT_PMOD;
  /* Some format has to be set anyway, and simple format would lose information */
  if(sim._type==DFT_FORMAT_PLAIN || sim._type==DFT_FORMAT_UNKNOWN)
    sim._type=DFT_FORMAT_STANDARD;
  /* Write file */
  if(dftWrite(&sim, sfile)) {
    fprintf(stderr, "Error in writing '%s': %s\n", sfile, dfterrmsg);
    dftEmpty(&sim); return(11);
  }
  if(verbose>0) fprintf(stdout, "simulated TAC(s) written in %s\n", sfile);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/