simoxygen.c File Reference

Simulation of [O-15]O2 tissue kinetics. More...

#include "tpcclibConfig.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include "tpccm.h"

Go to the source code of this file.

Functions
int	simOxygen (double *t, double *ca1, double *ca2, double *ca1i, double *ca2i, const int n, const double k1a, const double k2a, const double km, const double k1b, const double k2b, const double vb, const double fa, const int vvm, double *scpet, double *sct1, double *sct2, double *sctab, double *sctvb1, double *sctvb2, double *scvb1, double *scvb2, const int verbose)

Detailed Description

Simulation of [O-15]O2 tissue kinetics.

Definition in file simoxygen.c.

Function Documentation

simOxygen()

int simOxygen	(	double *	t,
		double *	ca1,
		double *	ca2,
		double *	ca1i,
		double *	ca2i,
		const int	n,
		const double	k1a,
		const double	k2a,
		const double	km,
		const double	k1b,
		const double	k2b,
		const double	vb,
		const double	fa,
		const int	vvm,
		double *	scpet,
		double *	sct1,
		double *	sct2,
		double *	sctab,
		double *	sctvb1,
		double *	sctvb2,
		double *	scvb1,
		double *	scvb2,
		const int	verbose )

Simulate tissue and venous blood TACs using dual-input compartment model for [O-15]O2 (one tissue compartment for [O-15]O2, and another tissue compartment for its metabolite [O-15]H2O). Based on Mintun et al. J Nucl Med. 1984;25(2):177-187.

The units of rate constants must be related to the time unit of the data; 1/min and min, or 1/sec and sec.

Returns

Function returns 0 when successful, else a value >= 1.

Author

Vesa Oikonen

See also

simMBF, simC1, simC2, simC1DI, simDispersion, simDelay, simSamples

Parameters

t	Array of sample times.
ca1	Array of arterial blood activities of tracer1 ([O-15]O2).
ca2	Array of arterial blood activities of tracer2 ([O-15]H2O).
ca1i	Array of AUC 0-t of arterial tracer1 activities; NULL if not available.
ca2i	Array of AUC 0-t of arterial tracer2 activities; NULL if not available.
n	Nr of samples (array lengths).
k1a	Rate constant of the model for tracer1 (from blood to C1).
k2a	Rate constant of the model for tracer1 (from C1 to blood).
km	Rate constant of the model (from tracer1 in C1 to tracer2 in C2).
k1b	Rate constant of the model for tracer2 (from blood to C2).
k2b	Rate constant of the model for tracer2 (from C2 to blood).
vb	Vascular volume fraction [0-1].
fa	Arterial fraction of vascular volume [0-1].
vvm	Vascular volume modelling: set to 0 to use Cpet = Vb*Cb + (1-Vb)*Ct, or set to 1 to use Cpet = Vb*Cb + Ct.
scpet	Pointer for TTAC array to be simulated; allocate in the calling program or set to NULL if not needed.
sct1	Simulated TAC of tracer1 in tissue; allocate in the calling program or set to NULL if not needed.
sct2	Simulated TAC of tracer2 in tissue; allocate in the calling program or set to NULL if not needed.
sctab	Total arterial contribution to PET TTAC; allocate in the calling program or set to NULL if not needed.
sctvb1	Venous tracer1 contribution to PET TAC; allocate in the calling program or set to NULL if not needed.
sctvb2	Venous tarcer1 contribution to PET TAC; allocate in the calling program or set to NULL if not needed.
scvb1	Venous BTAC of tracer1; allocate in the calling program or set to NULL if not needed.
scvb2	Venous BTAC of tracer2; allocate in the calling program or set to NULL if not needed.
verbose	Verbose level; if zero, then nothing is printed into stdout or stderr.

Definition at line 29 of file simoxygen.c.

  {
  if(verbose>0) {
    printf("%s()\n", __func__);
    if(verbose>1) {
      printf("  k1a := %g\n", k1a);
      printf("  k2a := %g\n", k2a);
      printf("  km := %g\n", km);
      printf("  k1b := %g\n", k1b);
      printf("  k2b := %g\n", k2b);
      printf("  vb := %g\n", vb);
      printf("  fa := %g\n", fa);
      printf("  n := %d\n", n);
    }
  }
 
  /* Check for data */
  if(n<2) return 1;
  if(t==NULL) return 2;
  if(ca1==NULL || ca2==NULL) return 3;
 
  /* Check parameters */
  if(k1a<0.0 || k1b<0.0 || k2a<0.0 || k2b<0.0) return(4);
  if(vb<0.0 || vb>=1.0) return(5);
  if(fa<0.0 || fa>1.0) return(6);
  double va=fa*vb;       // arterial volume fraction in tissue
  double vv=(1.0-fa)*vb; // venous volume fraction in tissue
 
  /* Set initial condition */
  double t_last=0.0; // previous sample time
  if(t[0]<t_last) t_last=t[0];
  /* Concentrations, integrals, and previous values are zero */
  double cba1i=0.0, cba1_last=0.0;
  double cba2i=0.0, cba2_last=0.0;
  double ct1=0.0, ct1_last=0.0, ct1i=0.0, ct1i_last=0.0;
  double ct2=0.0, ct2_last=0.0, ct2i=0.0, ct2i_last=0.0;
  double cvb1=0.0, cvb2=0.0;
 
  /* Calculate curves */
  double p, q;
  double dt2; // delta t / 2
  for(int i=0; i<n; i++) {
    dt2=0.5*(t[i]-t_last);
    if(dt2>0.0) {
      /* arterial integrals */
      if(ca1i!=NULL) cba1i=ca1i[i]; else cba1i+=(ca1[i]+cba1_last)*dt2;
      if(ca2i!=NULL) cba2i=ca2i[i]; else cba2i+=(ca2[i]+cba2_last)*dt2;
      /* partial results */
      p=ct1i_last+dt2*ct1_last;
      q=ct2i_last+dt2*ct2_last;
      /* 1st tissue compartment and its integral */
      ct1 = (k1a*cba1i - (k2a+km)*p) / (1.0 + dt2*(k2a+km));
      ct1i = ct1i_last + dt2*(ct1_last+ct1);
      /* 2nd tissue compartment and its integral */
      ct2 = (km*ct1i + k1b*cba2i - k2b*q) / (1.0 + dt2*k2b);
      ct2i = ct2i_last + dt2*(ct2_last+ct2);
    }
    /* Venous BTACs */
    if(k1a>0.0 && k2a>0.0) cvb1=ct1/(k1a/k2a);
    else if(k2a>0.0) cvb1=0.0; else cvb1=ca1[i];
    if(k1b>0.0 && k2b>0.0) cvb2=ct2/(k1b/k2b);
    else if(k2b>0.0) cvb2=0.0; else cvb2=ca2[i];
    /* copy values to argument arrays */
    if(scpet!=NULL) {
      if(vvm==0) scpet[i]= va*(ca1[i]+ca2[i]) + vv*(cvb1+cvb2) + (1.0-vb)*(ct1+ct2);
      else scpet[i]= va*(ca1[i]+ca2[i]) + vv*(cvb1+cvb2) + (ct1+ct2);
    }
    if(sct1!=NULL) {
      sct1[i]=ct1; if(vvm==0) sct1[i]*=(1.0-vb);
    }
    if(sct2!=NULL) {
      sct2[i]=ct2; if(vvm==0) sct2[i]*=(1.0-vb);
    }
    if(sctab!=NULL) {
      sctab[i]=va*(ca1[i]+ca2[i]); 
    }
    if(sctvb1!=NULL) {
      sctvb1[i]=vv*cvb1; 
    }
    if(sctvb2!=NULL) {
      sctvb2[i]=vv*cvb2; 
    }
    if(scvb1!=NULL) scvb1[i]=cvb1;
    if(scvb2!=NULL) scvb2[i]=cvb2;
    /* prepare for the next loop */
    t_last=t[i]; 
    cba1_last=ca1[i]; cba2_last=ca2[i];
    ct1_last=ct1; ct1i_last=ct1i;
    ct2_last=ct2; ct2i_last=ct2i;
  } // next sample
 
  if(verbose>2) {
    printf("AUC 0-%g:\n", t_last);
    printf(" cba1i := %g\n", cba1i);
    printf(" cba2i := %g\n", cba2i);
    printf(" ct1i := %g\n", ct1i_last);
    printf(" ct2i := %g\n", ct2i_last);
  }
 
  return 0;
}