simdicm.c File Reference

Simulation of dual-input compartmental models. 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	simC1DI (double *t, double *cba, double *cbb, const int nr, const double k1a, const double k1b, const double k2, double *ct)
int	simC3DIvs (double *t, double *ca1, double *ca2, double *cb, const int nr, const double k1, const double k2, const double k3, const double k4, const double k5, const double k6, const double k1b, const double k2b, const double f, const double vb, const double fa, const int vvm, double *scpet, double *sct1, double *sct2, double *sct3, double *sct1b, double *sctab, double *sctvb)
int	simC4DIvp (double *t, double *ca1, double *ca2, double *cb, const int nr, const double k1, const double k2, const double k3, const double k4, const double k5, const double k6, const double k7, const double km, double k1b, const double k2b, const double f, const double vb, double fa, const int vvm, double *scpet, double *sct1, double *sct2, double *sct3, double *sct1b, double *sctab, double *sctvb, const int verbose)
int	simC4DIvs (double *t, double *ca1, double *ca2, double *cb, const int nr, const double k1, const double k2, const double k3, const double k4, const double k5, const double k6, const double k7, const double km, double k1b, double k2b, double f, double vb, double fa, const int vvm, double *scpet, double *sct1, double *sct2, double *sct3, double *sct1b, double *sctab, double *sctvb, const int verbose)

Detailed Description

Simulation of dual-input compartmental models.

Definition in file simdicm.c.

Function Documentation

simC1DI()

int simC1DI	(	double *	t,
		double *	cba,
		double *	cbb,
		const int	nr,
		const double	k1a,
		const double	k1b,
		const double	k2,
		double *	ct )

Simulate tissue TAC using dual-input tissue compartment model with a single tissue compartment.

Simulates tissue TAC using dual-input tissue compartment model (1 compartment, common for both inputs) at input TAC times. Vascular volume is not accounted for. The units of rate constants must be related to the time unit; 1/min and min, or 1/sec and sec.

Returns

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

Author

Vesa Oikonen

See also

simC1, simDispersion, simOxygen, simC3vp, simC3vs, parExamplePerfectBolus

Todo

Finish coding and test.

Parameters

t	Array of time values.
cba	Array of concentrations of input A.
cbb	Array of concentrations of input B.
nr	Number of values in TACs.
k1a	Rate constant of the model for input A (from blood to C1).
k1b	Rate constant of the model for input B (from blood to C1).
k2	Rate constant of the model (from C1 to blood).
ct	Pointer for TAC array to be simulated; must be allocated.

Definition at line 30 of file simdicm.c.

  {
  /* Check for data */
  if(nr<2 || t==NULL || cba==NULL || cbb==NULL) return(1);
  if(ct==NULL) return(2);
 
  /* Check parameters */
  if(!(k1a>=0.0)) return(3);
  if(!(k1b>=0.0)) return(4);
  if(!(k2>=0.0)) return(5);
 
  /* Calculate curves */
  double t_last=0.0; if(t[0]<t_last) t_last=t[0]; 
  double cbai=0.0, cba_last=0.0, cbbi=0.0, cbb_last=0.0;
  double ct1=0.0, ct1i=0.0, ct1_last=0.0, ct1i_last=0.0;
  for(int i=0; i<nr; i++) {
    /* delta time / 2 */
    double dt2=0.5*(t[i]-t_last);
    /* calculate values */
    if(dt2<0.0) {
      return(6);
    } else if(dt2>0.0) {
      /* input integrals */
      cbai+=(cba[i]+cba_last)*dt2;
      cbbi+=(cbb[i]+cbb_last)*dt2;
      /* tissue compartment and its integral */
      ct1 = (k1a*cbai + k1b*cbbi - k2*(ct1i_last+dt2*ct1_last)) / (1.0 + dt2*k2);
      ct1i = ct1i_last + dt2*(ct1_last+ct1);
    }
    /* copy values to argument array */
    ct[i]=ct1;
    /* prepare to the next loop */
    t_last=t[i]; cba_last=cba[i]; cbb_last=cbb[i];
    ct1_last=ct1; ct1i_last=ct1i;
  }
 
  return(0);
}

simC3DIvs()

int simC3DIvs	(	double *	t,
		double *	ca1,
		double *	ca2,
		double *	cb,
		const int	nr,
		const double	k1,
		const double	k2,
		const double	k3,
		const double	k4,
		const double	k5,
		const double	k6,
		const double	k1b,
		const double	k2b,
		const double	f,
		const double	vb,
		const double	fa,
		const int	vvm,
		double *	scpet,
		double *	sct1,
		double *	sct2,
		double *	sct3,
		double *	sct1b,
		double *	sctab,
		double *	sctvb )

Simulate tissue TAC using dual-input tissue compartment model with compartments in series.

Simulates tissue TAC using dual-input tissue compartment model (1-3 compartments in series for tracer1, and 1 compartment for tracer2) at plasma TAC times, considering also contribution of arterial and venous vasculature, but no exchange between compartments for tracer1 and tracer2. The units of rate constants must be related to the time unit; 1/min and min, or 1/sec and sec. If blood flow is set to 0, function assumes that f>>k1, and Cvb=Cab.

Returns

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

Author

Vesa Oikonen

See also

simC4DIvp, simC4DIvs, simC1DI, simC3vp, simC3vs, simDelay, parExamplePerfectBolus

Parameters

t	Array of time values.
ca1	Array of arterial plasma activities of tracer1.
ca2	Array of arterial plasma activities of tracer2.
cb	Array of arterial blood activities.
nr	Number of values in TACs.
k1	Rate constant of the model for tracer1 (from plasma to C1).
k2	Rate constant of the model for tracer1 (from C1 to plasma).
k3	Rate constant of the model for tracer1 (from C1 to C2).
k4	Rate constant of the model for tracer1 (from C2 to C1).
k5	Rate constant of the model for tracer1 (from C2 to C3).
k6	Rate constant of the model for tracer1 (from C3 to C2).
k1b	Rate constant of the model for tracer2 (from plasma to C4).
k2b	Rate constant of the model for tracer2 (from C4 to plasma).
f	Blood flow; if 0, function assumes that f>>k1, and Cvb=Cab.
vb	Vascular volume fraction.
fa	Arterial fraction of vascular volume.
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 TAC array to be simulated; must be allocated.
sct1	Pointer for 1st tracer1 compartment TAC, or NULL if not needed.
sct2	Pointer for 2nd tracer1 compartment TAC, or NULL if not needed.
sct3	Pointer for 3rd tracer1 compartment TAC, or NULL if not needed.
sct1b	Pointer for 1st tracer2 compartment TAC, or NULL if not needed.
sctab	Pointer for arterial TAC in tissue, or NULL if not needed.
sctvb	Pointer for venous TAC in tissue, or NULL if not needed.

Definition at line 101 of file simdicm.c.

  {
  int i;
  double b, c, d, e, w, z, dt2, va, vv;
  double ca1i, ca1_last, ca2i, ca2_last, t_last, dct, cvb;
  double ct1, ct1_last, ct2, ct2_last, ct3, ct3_last;
  double ct1i, ct1i_last, ct2i, ct2i_last, ct3i, ct3i_last;
  double ct1b, ct1b_last, ct1bi, ct1bi_last;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(scpet==NULL) return 2;
 
  /* Check parameters */
  if(k1<0.0) return 3;
  if(vb<0.0 || vb>=1.0) return 4;
  if(fa<=0.0 || fa>1.0) return 5;
  va=fa*vb; vv=(1.0-fa)*vb;
 
  /* Calculate curves */
  t_last=0.0; if(t[0]<t_last) t_last=t[0];
  ca1i=ca1_last=ca2i=ca2_last=0.0;
  ct1_last=ct2_last=ct3_last=ct1i_last=ct2i_last=ct3i_last=0.0;
  ct1b_last=ct1bi_last=0.0;
  ct1=ct2=ct3=ct1b=ct1i=ct2i=ct3i=ct1bi=0.0;
  for(i=0; i<nr; i++) {
    /* delta time / 2 */
    dt2=0.5*(t[i]-t_last);
    /* calculate values */
    if(dt2<0.0) {
      return 5;
    } else if(dt2>0.0) {
      /* arterial integrals */
      ca1i+=(ca1[i]+ca1_last)*dt2;
      ca2i+=(ca2[i]+ca2_last)*dt2;
      /* partial results */
      b=ct1i_last+dt2*ct1_last;
      c=ct2i_last+dt2*ct2_last;
      d=ct3i_last+dt2*ct3_last;
      e=ct1bi_last+dt2*ct1b_last;
      w=k4 + k5 - (k5*k6*dt2)/(1.0+k6*dt2);
      z=1.0+w*dt2;
      /* 1st tissue compartment and its integral */
      ct1 = (
          + k1*z*ca1i + (k3*k4*dt2 - (k2+k3)*z)*b
          + k4*c + k4*k6*dt2*d/(1.0+k6*dt2)
        ) / ( z*(1.0 + dt2*(k2+k3)) - k3*k4*dt2*dt2 );
      ct1i = ct1i_last + dt2*(ct1_last+ct1);
      ct1b = (k1b*ca2i - k2b*e) / (1.0 + dt2*k2b);
      ct1bi = ct1bi_last + dt2*(ct1b_last+ct1b);
      /* 2nd tissue compartment and its integral */
      ct2 = (k3*ct1i - w*c + k6*d/(1.0+k6*dt2)) / z;
      ct2i = ct2i_last + dt2*(ct2_last+ct2);
      /* 3rd tissue compartment and its integral */
      ct3 = (k5*ct2i - k6*d) / (1.0 + k6*dt2);
      ct3i = ct3i_last + dt2*(ct3_last+ct3);
    }
    /* Venous curve */
    if(f>0.) {
      dct = k1*ca1[i] - k2*ct1 + k1b*ca2[i] - k2b*ct1b;
      cvb = cb[i] - dct/f;
    } else cvb=cb[i];
    /* copy values to argument arrays */
    if(vvm==0) scpet[i]= va*cb[i] + vv*cvb + (1.0-vb)*(ct1+ct2+ct3+ct1b);
    else scpet[i]= va*cb[i] + vv*cvb + (ct1+ct2+ct3+ct1b);
    if(sct1!=NULL) {sct1[i]=(1.0-vb)*ct1; if(vvm==0) sct1[i]*=(1.0-vb);}
    if(sct2!=NULL) {sct2[i]=(1.0-vb)*ct2; if(vvm==0) sct2[i]*=(1.0-vb);}
    if(sct3!=NULL) {sct3[i]=(1.0-vb)*ct3; if(vvm==0) sct3[i]*=(1.0-vb);}
    if(sct1b!=NULL) {sct1b[i]=(1.0-vb)*ct1b; if(vvm==0) sct1b[i]*=(1.0-vb);}
    if(sctab!=NULL) sctab[i]=va*cb[i];
    if(sctvb!=NULL) sctvb[i]=vv*cvb;
    /* prepare to the next loop */
    t_last=t[i]; ca1_last=ca1[i]; ca2_last=ca2[i];
    ct1_last=ct1; ct1i_last=ct1i;
    ct2_last=ct2; ct2i_last=ct2i;
    ct3_last=ct3; ct3i_last=ct3i;
    ct1b_last=ct1b; ct1bi_last=ct1bi;
  }
 
  return 0;
}

simC4DIvp()

int simC4DIvp	(	double *	t,
		double *	ca1,
		double *	ca2,
		double *	cb,
		const int	nr,
		const double	k1,
		const double	k2,
		const double	k3,
		const double	k4,
		const double	k5,
		const double	k6,
		const double	k7,
		const double	km,
		double	k1b,
		const double	k2b,
		const double	f,
		const double	vb,
		double	fa,
		const int	vvm,
		double *	scpet,
		double *	sct1,
		double *	sct2,
		double *	sct3,
		double *	sct1b,
		double *	sctab,
		double *	sctvb,
		const int	verbose )

Simulate tissue TAC using dual-input tissue compartment model (compartments 2 and 3 in parallel for tracer1, and 1 compartment for tracer2) at plasma TAC sample times, considering also contribution of arterial and venous vasculature, and transfer of tracer1 to tracer2.

The units of rate constants must be related to the time unit; 1/min and min, or 1/sec and sec. If blood flow is set to 0, function assumes that f>>k1, and Cvb=Cab. Reference: TPCMOD0001 Appendix C. Tested with program p2t_di -parallel.

Returns

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

Author

Vesa Oikonen

See also

simC3DIvs, simC4DIvs, simC1DI, simC3vp, simC3vs, parExamplePerfectBolus

Parameters

t	Array of time values.
ca1	Array of arterial plasma activities of tracer1 (parent tracer).
ca2	Array of arterial plasma activities of tracer2 (metabolite).
cb	Array of (total) arterial blood activities.
nr	Number of values in TACs.
k1	Rate constant of the model for tracer1 (from plasma to C1).
k2	Rate constant of the model for tracer1 (from C1 to plasma).
k3	Rate constant of the model for tracer1 (from C1 to C2).
k4	Rate constant of the model for tracer1 (from C2 to C1).
k5	Rate constant of the model for tracer1 (from C1 to C3).
k6	Rate constant of the model for tracer1 (from C3 to C1).
k7	Rate constant of the model for tracer1 (from C3 to plasma).
km	Rate constant of the model (from tracer1 in C1 to tracer2 in C4).
k1b	Rate constant of the model for tracer2 (from plasma to C4).
k2b	Rate constant of the model for tracer2 (from C4 to plasma).
f	Blood flow; if 0, function assumes that f>>k1, and Cvb=Cab.
vb	Vascular volume fraction (0<=Vb<1).
fa	Arterial fraction of vascular volume (0<=fa<=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 TAC array to be simulated; must be allocated.
sct1	Pointer for 1st tracer1 compartment TAC, or NULL if not needed.
sct2	Pointer for 2nd tracer1 compartment TAC, or NULL if not needed.
sct3	Pointer for 3rd tracer1 compartment TAC, or NULL if not needed.
sct1b	Pointer for 1st tracer2 compartment TAC, or NULL if not needed.
sctab	Pointer for arterial TAC in tissue, or NULL if not needed.
sctvb	Pointer for venous TAC in tissue, or NULL if not needed.
verbose	Verbose level; if zero, then nothing is printed into stdout or stderr.

Definition at line 251 of file simdicm.c.

  {
  int i;
  double b, c, d, e, pt, qt, dt2, va, vv;
  double ca1i, ca1_last, ca2i, ca2_last, t_last, dct, cvb;
  double ct1, ct1_last, ct2, ct2_last, ct3, ct3_last;
  double ct1i, ct1i_last, ct2i, ct2i_last, ct3i, ct3i_last;
  double ct1b, ct1b_last, ct1bi, ct1bi_last;
 
 
  if(verbose>0) {
    printf("%s()\n", __func__);
    if(verbose>1) {
      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("  k7 := %g\n", k7);
      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("  f := %g\n", f);
    }
  }
 
  /* Check for data */
  if(nr<2) return 1;
  if(scpet==NULL) return 2;
 
  /* Check parameters */
  if(k1<0.0 || k1b<0.0) return 3;
  if(vb<0.0 || vb>=1.0) return 4;
  if(fa<0.0 || fa>1.0) return 5;
  va=fa*vb; vv=(1.0-fa)*vb;
 
  /* Calculate curves */
  t_last=0.0; if(t[0]<t_last) t_last=t[0];
  ca1i=ca1_last=ca2i=ca2_last=0.0;
  ct1_last=ct2_last=ct3_last=ct1i_last=ct2i_last=ct3i_last=ct1b_last=
  ct1bi_last=0.0;
  ct1=ct2=ct3=ct1b=ct1i=ct2i=ct3i=ct1bi=0.0;
  for(i=0; i<nr; i++) {
    /* delta time / 2 */
    dt2=0.5*(t[i]-t_last);
    /* calculate values */
    if(dt2<0.0) {
      return 5;
    } else if(dt2>0.0) {
      /* arterial integrals */
      ca1i+=(ca1[i]+ca1_last)*dt2;
      ca2i+=(ca2[i]+ca2_last)*dt2;
      /* partial results */
      b=ct1i_last+dt2*ct1_last;
      c=ct2i_last+dt2*ct2_last;
      d=ct3i_last+dt2*ct3_last;
      e=ct1bi_last+dt2*ct1b_last;
      pt=k6+k7;
      qt=k2+k3+k5+km-(k3*k4*dt2)/(1.0+k4*dt2)-(k5*k6*dt2)/(1.0+pt*dt2);
      /* 1st tissue compartment and its integral */
      ct1 = (k1/(1.0+qt*dt2))*ca1i
          - (qt/(1.0+qt*dt2))*b
          + (k4/((1.0+qt*dt2)*(1.0+k4*dt2)))*c
          + (k6/((1.0+qt*dt2)*(1.0+pt*dt2)))*d;
      ct1i = ct1i_last + dt2*(ct1_last+ct1);
      /* 2nd tissue compartment and its integral */
      ct2 = (k3/(1.0+k4*dt2))*ct1i
          - (k4/(1.0+k4*dt2))*c;
      ct2i = ct2i_last + dt2*(ct2_last+ct2);
      /* 3rd tissue compartment and its integral */
      ct3 = (k5/(1.0+pt*dt2))*ct1i
          - (pt/(1.0+pt*dt2))*d;
      ct3i = ct3i_last + dt2*(ct3_last+ct3);
      /* 4th tissue compartment (the 1st for tracer 2) and its integral */
      ct1b = (k1b/(1.0+k2b*dt2))*ca2i
           - (k2b/(1.0+k2b*dt2))*e
           + (km/(1.0+k2b*dt2))*ct1i;
      ct1bi = ct1bi_last + dt2*(ct1b_last+ct1b);
    }
    /* Venous curve */
    if(f>0.) {
      dct = k1*ca1[i] - k2*ct1 - k7*ct3 + k1b*ca2[i] - k2b*ct1b;
      cvb = cb[i] - dct/f;
    } else cvb=cb[i];
    /* copy values to argument arrays */
    if(vvm==0) scpet[i]= va*cb[i] + vv*cvb + (1.0-vb)*(ct1+ct2+ct3+ct1b);
    else scpet[i]= va*cb[i] + vv*cvb + (ct1+ct2+ct3+ct1b);
    if(sct1!=NULL) {sct1[i]=(1.0-vb)*ct1; if(vvm==0) sct1[i]*=(1.0-vb);}
    if(sct2!=NULL) {sct2[i]=(1.0-vb)*ct2; if(vvm==0) sct2[i]*=(1.0-vb);}
    if(sct3!=NULL) {sct3[i]=(1.0-vb)*ct3; if(vvm==0) sct3[i]*=(1.0-vb);}
    if(sct1b!=NULL) {sct1b[i]=(1.0-vb)*ct1b; if(vvm==0) sct1b[i]*=(1.0-vb);}
    if(sctab!=NULL) sctab[i]=va*cb[i];
    if(sctvb!=NULL) sctvb[i]=vv*cvb;
    /* prepare to the next loop */
    t_last=t[i]; ca1_last=ca1[i]; ca2_last=ca2[i];
    ct1_last=ct1; ct1i_last=ct1i;
    ct2_last=ct2; ct2i_last=ct2i;
    ct3_last=ct3; ct3i_last=ct3i;
    ct1b_last=ct1b; ct1bi_last=ct1bi;
  }
  
  if(verbose>2) {
    printf("AUC 0-%g:\n", t_last);
    printf(" ca1i := %g\n", ca1i);
    printf(" ca2i := %g\n", ca2i);
    printf(" ct1i := %g\n", ct1i_last);
    printf(" ct2i := %g\n", ct2i_last);
    printf(" ct3i := %g\n", ct3i_last);
    printf(" ct1bi := %g\n", ct1bi_last);
  }
 
  return 0;
}

simC4DIvs()

int simC4DIvs	(	double *	t,
		double *	ca1,
		double *	ca2,
		double *	cb,
		const int	nr,
		const double	k1,
		const double	k2,
		const double	k3,
		const double	k4,
		const double	k5,
		const double	k6,
		const double	k7,
		const double	km,
		double	k1b,
		double	k2b,
		double	f,
		double	vb,
		double	fa,
		const int	vvm,
		double *	scpet,
		double *	sct1,
		double *	sct2,
		double *	sct3,
		double *	sct1b,
		double *	sctab,
		double *	sctvb,
		const int	verbose )

Simulate tissue TAC using dual-input tissue compartment model (1-3 compartments in series for tracer1, and 1 compartment for tracer2) at plasma TAC times, considering also contribution of arterial and venous vasculature, and transfer of tracer1 to tracer2.

The units of rate constants must be related to the time unit; 1/min and min, or 1/sec and sec. If blood flow is set to 0, function assumes that f>>k1, and Cvb=Cab. Reference: TPCMOD0001 Appendix B. Tested with program p2t_di -series.

Returns

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

Author

Vesa Oikonen

See also

simC3DIvs, simC4DIvp, simC1DI, simC3vp, simC3vs, parExamplePerfectBolus

Parameters

t	Array of time values.
ca1	Array of arterial plasma activities of tracer1 (parent tracer).
ca2	Array of arterial plasma activities of tracer2 (metabolite).
cb	Array of (total) arterial blood activities.
nr	Number of values in TACs.
k1	Rate constant of the model for tracer1 (from plasma to C1).
k2	Rate constant of the model for tracer1 (from C1 to plasma).
k3	Rate constant of the model for tracer1 (from C1 to C2).
k4	Rate constant of the model for tracer1 (from C2 to C1).
k5	Rate constant of the model for tracer1 (from C2 to C3).
k6	Rate constant of the model for tracer1 (from C3 to C2).
k7	Rate constant of the model for tracer1 (from C3 to plasma).
km	Rate constant of the model (from tracer1 in C1 to tracer2 in C4).
k1b	Rate constant of the model for tracer2 (from plasma to C4).
k2b	Rate constant of the model for tracer2 (from C4 to plasma).
f	Blood flow; if 0, function assumes that f>>k1, and Cvb=Cab.
vb	Vascular volume fraction.
fa	Arterial fraction of vascular volume.
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 TAC array to be simulated; must be allocated.
sct1	Pointer for 1st tracer1 compartment TAC, or NULL if not needed.
sct2	Pointer for 2nd tracer1 compartment TAC, or NULL if not needed.
sct3	Pointer for 3rd tracer1 compartment TAC, or NULL if not needed.
sct1b	Pointer for 1st tracer2 compartment TAC, or NULL if not needed.
sctab	Pointer for arterial TAC in tissue, or NULL if not needed.
sctvb	Pointer for venous TAC in tissue, or NULL if not needed.
verbose	Verbose level; if zero, then nothing is printed into stdout or stderr.

Definition at line 441 of file simdicm.c.

  {
  int i;
  double b, c, d, e, pt, qt, rt, dt2, va, vv;
  double ca1i, ca1_last, ca2i, ca2_last, t_last, dct, cvb;
  double ct1, ct1_last, ct2, ct2_last, ct3, ct3_last;
  double ct1i, ct1i_last, ct2i, ct2i_last, ct3i, ct3i_last;
  double ct1b, ct1b_last, ct1bi, ct1bi_last;
 
 
  if(verbose>0) {
    printf("%s()\n", __func__);
    if(verbose>1) {
      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("  k7 := %g\n", k7);
      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("  f := %g\n", f);
    }
  }
 
  /* Check for data */
  if(nr<2) return 1;
  if(scpet==NULL) return 2;
 
  /* Check parameters */
  if(k1<0.0 || k1b<0.0) return 3;
  if(vb<0.0 || vb>=1.0) return 4;
  if(fa<0.0 || fa>1.0) return 5;
  va=fa*vb; vv=(1.0-fa)*vb;
 
  /* Calculate curves */
  t_last=0.0; if(t[0]<t_last) t_last=t[0];
  ca1i=ca1_last=ca2i=ca2_last=0.0;
  ct1_last=ct2_last=ct3_last=ct1i_last=ct2i_last=ct3i_last=ct1b_last=
  ct1bi_last=0.0;
  ct1=ct2=ct3=ct1b=ct1i=ct2i=ct3i=ct1bi=0.0;
  for(i=0; i<nr; i++) {
    /* delta time / 2 */
    dt2=0.5*(t[i]-t_last);
    /* calculate values */
    if(dt2<0.0) {
      return 5;
    } else if(dt2>0.0) {
      /* arterial integrals */
      ca1i+=(ca1[i]+ca1_last)*dt2;
      ca2i+=(ca2[i]+ca2_last)*dt2;
      //printf("ca1[%d]=%g int=%g\n", i, ca1[i], ca1i);
      //printf("ca2[%d]=%g int=%g\n", i, ca2[i], ca2i);
      /* partial results */
      b=ct1i_last+dt2*ct1_last;
      c=ct2i_last+dt2*ct2_last;
      d=ct3i_last+dt2*ct3_last;
      e=ct1bi_last+dt2*ct1b_last;
      pt=k6+k7;
      qt=k4+k5-(k5*k6*dt2)/(1.0+pt*dt2);
      rt=k2+k3+km-(k3*k4*dt2)/(1.0+qt*dt2);
      /* 1st tissue compartment and its integral */
      ct1 = (k1/(1.0+rt*dt2))*ca1i
          - (rt/(1.0+rt*dt2))*b
          + (k4/((1.0+qt*dt2)*(1.0+rt*dt2)))*c
          + ((k4*k6*dt2)/((1.0+pt*dt2)*(1.0+qt*dt2)*(1.0+rt*dt2)))*d;
      ct1i = ct1i_last + dt2*(ct1_last+ct1);
      /* 2nd tissue compartment and its integral */
      ct2 = (k3/(1.0+qt*dt2))*ct1i
          - (qt/(1.0+qt*dt2))*c
          + (k6/((1.0+pt*dt2)*(1.0+qt*dt2)))*d;
      ct2i = ct2i_last + dt2*(ct2_last+ct2);
      /* 3rd tissue compartment and its integral */
      ct3 = (k5/(1.0+pt*dt2))*ct2i
          - (pt/(1.0+pt*dt2))*d;
      ct3i = ct3i_last + dt2*(ct3_last+ct3);
      /* 4th tissue compartment (the 1st for tracer 2) and its integral */
      ct1b = (k1b/(1.0+k2b*dt2))*ca2i
           - (k2b/(1.0+k2b*dt2))*e
           + (km/(1.0+k2b*dt2))*ct1i;
      ct1bi = ct1bi_last + dt2*(ct1b_last+ct1b);
    }
    /* Venous curve */
    if(f>0.) {
      dct = k1*ca1[i] - k2*ct1 - k7*ct3 + k1b*ca2[i] - k2b*ct1b;
      cvb = cb[i] - dct/f;
    } else cvb=cb[i];
    /* copy values to argument arrays; set very small values to zero */
    if(vvm==0) scpet[i]= va*cb[i] + vv*cvb + (1.0-vb)*(ct1+ct2+ct3+ct1b);
    else scpet[i]= va*cb[i] + vv*cvb + (ct1+ct2+ct3+ct1b);
    if(sct1!=NULL) {sct1[i]=(1.0-vb)*ct1; if(vvm==0) sct1[i]*=(1.0-vb);}
    if(sct2!=NULL) {sct2[i]=(1.0-vb)*ct2; if(vvm==0) sct2[i]*=(1.0-vb);}
    if(sct3!=NULL) {sct3[i]=(1.0-vb)*ct3; if(vvm==0) sct3[i]*=(1.0-vb);}
    if(sct1b!=NULL) {sct1b[i]=(1.0-vb)*ct1b; if(vvm==0) sct1b[i]*=(1.0-vb);}
    if(sctab!=NULL) sctab[i]=va*cb[i];
    if(sctvb!=NULL) sctvb[i]=vv*cvb;
    /* prepare to the next loop */
    t_last=t[i]; ca1_last=ca1[i]; ca2_last=ca2[i];
    ct1_last=ct1; ct1i_last=ct1i;
    ct2_last=ct2; ct2i_last=ct2i;
    ct3_last=ct3; ct3i_last=ct3i;
    ct1b_last=ct1b; ct1bi_last=ct1bi;
  }
  
  if(verbose>2) {
    printf("AUC 0-%g:\n", t_last);
    printf(" ca1i := %g\n", ca1i);
    printf(" ca2i := %g\n", ca2i);
    printf(" ct1i := %g\n", ct1i_last);
    printf(" ct2i := %g\n", ct2i_last);
    printf(" ct3i := %g\n", ct3i_last);
    printf(" ct1bi := %g\n", ct1bi_last);
  }
 
  return 0;
}