p2t_v3c.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[] = {
  "Deprecated, use sim_3tcm instead.",
  " ",
  "Simulation of PET tissue time-radioactivity concentration curve (TAC)",
  "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] plasmafile bloodfile K1 k2 k3 k4 k5 k6 f Vb fA 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, default)",
  "     or not written (-nosub) into the simfile.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Enter blood flow (f) in units (ml/(min*ml)); set f=0 to assume that f>>K1",
  "and Cvb=Cab. Enter the vascular volume fraction (Vb) and its arterial",
  "fraction (fA) as percentages.",
  "If the times in input files 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, k5 or k3 can be set to 0.",
  " ",
  "Example:",
  "     @P -nosub plasma.dat blood.dat 0.3 0.2 0.1 0.2 0 0 0.6 5 30 sim.dat",
  " ",
  "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",
  "  4) Activity concentration in 2nd tissue compartment, (1-Vb)*C2",
  "  5) Activity concentration in 3rd tissue compartment, (1-Vb)*C3",
  "  6) Arterial contribution to tissue activity, Vb*fA*Cab",
  "  7) Venous contribution to tissue activity, Vb*(1-fA)*Cvb",
  " ",
  "References:",
  "1. TPCMOD0001.",
  " ",
  "See also: sim_3tcm, p2t_3c, p2t_di, tacadd, tacren, simframe, 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;
  int        ret;
  int        fi, ri;
  int        save_only_total=0;
  int        parallel=0;
  DFT        plasma, sim;
  double     k1, k2, k3, k4, k5, k6, f, Vb, fA;
  double    *t, *ca, *cab, *cpet, *ct1, *ct2, *ct3, *ctab, *ctvb;
  char       pfile[FILENAME_MAX], bfile[FILENAME_MAX], sfile[FILENAME_MAX];
  char       tmp[FILENAME_MAX+128], *cptr;
 
 
  
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  dftInit(&plasma); dftInit(&sim);
  pfile[0]=bfile[0]=sfile[0]=(char)0;
  k1=k2=k3=k4=k5=k6=f=Vb=fA=-1.0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    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;
    }
    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(!bfile[0]) {
      strcpy(bfile, argv[ai]); continue;
    } else if(k1<0.0) {
      k1=atof_dpi(argv[ai]); if(k1>0.0) continue;
    } else if(k2<0.0) {
      k2=atof_dpi(argv[ai]); if(k2>=0.0) continue;
    } else if(k3<0.0) {
      k3=atof_dpi(argv[ai]); if(k3>=0.0) continue;
    } else if(k4<0.0) {
      k4=atof_dpi(argv[ai]); if(k4>=0.0) continue;
    } else if(k5<0.0) {
      k5=atof_dpi(argv[ai]); if(k5>=0.0) continue;
    } else if(k6<0.0) {
      k6=atof_dpi(argv[ai]); if(k6>=0.0) continue;
    } else if(f<0.0) {
      f=atof_dpi(argv[ai]);  if(f>=0.0)  continue;
    } else if(Vb<0.0) {
      Vb=0.01*atof_dpi(argv[ai]); if(Vb<0.0) Vb=0.0; continue;
    } else if(fA<0.0) {
      fA=0.01*atof_dpi(argv[ai]); if(fA<1.0E-012) fA=0.0; 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]) {tpcPrintUsage(argv[0], info, stdout); return(1);}
  if(k1<=0.0) {fprintf(stderr, "Error: k1 must be > 0.\n"); return(1);}
  if(f>0.0 && f<k1) {
    fprintf(stderr, "Error: f cannot be lower than K1!\n"); return(1);}
  if(Vb>=1.0 || fA<=0.0 || fA>1.0) {
    fprintf(stderr, "Error: invalid value of Vb!\n"); return(1);}
  if(fA<=0.0 || fA>1.0) {
    fprintf(stderr, "Error: invalid value of fA!\n"); return(1);}
  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) {
    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("bfile := %s\n", bfile);
    printf("sfile := %s\n", sfile);
    printf("parameters := %g %g %g %g %g %g\n", k1, k2, k3, k4, k5, k6);
    printf("f  := %g\n", f);
    printf("Vb := %g\n", Vb);
    printf("fA := %g\n", fA);
    printf("save_only_total := %d\n", save_only_total);
    printf("parallel := %d\n", parallel);
  }
 
 
  /*
   *  Read plasma data
   */
  if(verbose>1) printf("reading plasma\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");
  }
 
  /*
   *  Read blood data; interpolate it to the plasma times
   */
  if(verbose) 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[1].voiname, "Blood");
  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,
      NULL, NULL, plasma.frameNr);
  dftEmpty(&sim);
  if(ret) {
    fprintf(stderr, "Error: cannot interpolate blood data.\n");
    dftEmpty(&plasma); return(3);
  }
  if(verbose>9) dftPrint(&plasma);
 
 
  /*
   *  Allocate memory for simulated tissue data.
   */
  if(verbose>1) printf("allocating memory\n");
  /* allocate */
  ret=dftSetmem(&sim, plasma.frameNr, 6);
  if(ret) {
    fprintf(stderr, 
            "Error (%d): cannot allocate memory for simulated TACs.\n", ret);
    dftEmpty(&plasma); return(4);
  }
  sim.frameNr=plasma.frameNr; sim.voiNr=6;
  /* 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, "Cab");
        sim.voi[ri].size=100.0*fA*Vb;
        break;
      case 5: 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; ca=plasma.voi[0].y; cab=plasma.voi[1].y;
  cpet=sim.voi[0].y; ct1=sim.voi[1].y; ct2=sim.voi[2].y; ct3=sim.voi[3].y;
  ctab=sim.voi[4].y; ctvb=sim.voi[5].y;
  if(parallel==0)
    ret=simC3vs(t, ca, cab, sim.frameNr, k1, k2, k3, k4, k5, k6, f, Vb, fA,
                cpet, ct1, ct2, ct3, ctab, ctvb);
  else
    ret=simC3vp(t, ca, cab, sim.frameNr, k1, k2, k3, k4, k5, k6, f, Vb, fA,
                cpet, ct1, ct2, ct3, ctab, ctvb);
  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;
  /* Set comments */
  dftSetComments(&sim);
  sprintf(tmp, "# plasmafile := %s\n", pfile); strcat(sim.comments, tmp);
  sprintf(tmp, "# bloodfile := %s\n", bfile); strcat(sim.comments, tmp);
  strcpy(tmp, "# model := ");
  if(parallel==0) strcat(tmp, "C3VS\n"); else strcat(tmp, "C3VP\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, "# 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);
  dftEmpty(&sim);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/