simoxygen.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
/*****************************************************************************/
#include "tpccm.h"
/*****************************************************************************/
 
/*****************************************************************************/
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
) {
  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;
}
/*****************************************************************************/
 
/*****************************************************************************/