tpccm.h File Reference

Header file for libtpccm. More...

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

Go to the source code of this file.

Data Structures
struct	ICMPARC

Functions
int	simMBF (double *t, double *ci, const int nr, const double k1, const double k2, const double Vfit, double *ct)
int	simC1 (double *t, double *ca, const int nr, const double k1, const double k2, double *ct)
int	simC1_i (double *t, double *cai, const int nr, const double k1, const double k2, double *ct)
int	simC1_d (double *t, double *ca, const int nr, const double k1, const double k2, double *ct)
int	simC2 (double *t, double *ca, const int nr, const double k1, const double k2, const double k3, const double k4, double *ct, double *cta, double *ctb)
int	simC2_i (double *t, double *cai, const int nr, const double k1, const double k2, const double k3, const double k4, double *ct, double *cta, double *ctb)
int	simC3s (double *t, double *ca, const int nr, double k1, double k2, double k3, double k4, double k5, double k6, double *ct, double *cta, double *ctb, double *ctc)
int	simC3vs (double *t, double *ca, 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 f, const double vb, const double fa, const int vvm, double *cpet, double *cta, double *ctb, double *ctc, double *ctab, double *ctvb)
int	simC3p (double *t, double *ca, const int nr, const double k1, const double k2, const double k3, const double k4, const double k5, const double k6, double *ct, double *cta, double *ctb, double *ctc)
int	simC3vp (double *t, double *ca, 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 f, const double vb, const double fa, const int vvm, double *cpet, double *cta, double *ctb, double *ctc, double *ctab, double *ctvb)
int	simC2l (double *t, double *ca, const int nr, const double k1, const double k2, const double k3, const double kLoss, double *ct, double *cta, double *ctb)
int	simC2vl (double *t, double *ca, double *cb, const int nr, const double k1, const double k2, const double k3, const double kL, const double f, const double vb, const double fa, const int vvm, double *cpet, double *cta, double *ctb, double *ctab, double *ctvb)
int	simC3vpKLoss (double *t, double *ca, 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 kLoss, const double f, const double vb, const double fa, const int vvm, double *cpet, double *cta, double *ctb, double *ctc, double *ctab, double *ctvb)
int	simRTCM (double *t, double *cr, const int nr, const double R1, const double k2, const double k3, const double k4, double *ct, double *cta, double *ctb)
int	simSRTM (double *t, double *cr, const int nr, const double R1, const double k2, const double BP, double *ct)
int	simTRTM (double *t, double *cr, const int nr, const double R1, const double k2, double k3, double *ct)
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, 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, 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, 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, const int verbose)
int	simDispersion (double *x, double *y, const int n, const double tau1, const double tau2, double *tmp)
int	corDispersion (double *x, double *y, const int n, const double tau, double *tmp)
int	simTTM (double *t, double *c0, const int n, const double k, const int cn, double *cout)
int	simTTM_i (double *t, double *c0i, const int n, const double k, const int cn, double *cout)
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)
int	convolve1D (double *data, const int n, double *kernel, const int m, double *out)
	Calculates the convolution sum of a discrete real data set data[0..n-1] and a discretized response function kernel[0..m].
int	simIsSteadyInterval (double *x, const int n, double *f)
void	icmparcInit (ICMPARC *d)
int	icmparcAddMetabolites (ICMPARC *d, unsigned const int mNr)
int	icmparcAllocateTACs (ICMPARC *d, unsigned const int sNr, const int sub)
void	icmparcFree (ICMPARC *d)
int	simBTAC (double *t, const unsigned int nr, ICMPARC *p, double *cb)
int	simBM1 (double *t, double *ct, const int nr, const double k, double *cp)

Detailed Description

Header file for libtpccm.

Header file for compartmental model library.

Author

Vesa Oikonen

Copyright

(c) Turku PET Centre

Definition in file tpccm.h.

Function Documentation

convolve1D()

int convolve1D	(	double *	data,
		const int	n,
		double *	kernel,
		const int	m,
		double *	out )

Calculates the convolution sum of a discrete real data set data[0..n-1] and a discretized response function kernel[0..m].

Convolution is not aware of the step size (default is 1); if step size is not 1 (it usually isn't), the step size must be taken into account either when computing the kernel or by scaling the convolution sum.

Remarks

This function does not use FFT, which would be faster with very large dataset, but slightly less precise.

See also

tacInterpolateToEqualLengthFrames, simDispersion, simC1

Returns

Function returns 0 when successful, or 1 if input data is not valid.

Parameters

data	Data array of length n-1 to be convolved, including any user-defined zero-padding.
n	Nr of data values.
kernel	Response function values in an array of length m.
m	Length of kernel array.
out	The convolved sum data is returned in out[0..n-1]; this must not overlap the input data.

Definition at line 27 of file convolut.c.

  {
  if(n<1 || m<1 || data==NULL || kernel==NULL || out==NULL || out==data) return(1);
 
  for(int di=0; di<n; di++) out[di]=0.0;
 
  for(int di=m-1; di<n; di++) {
    int dj=di;
    for(int k=0; k<m; k++)
      out[di] += data[dj--] * kernel[k];
  }
 
  for(int di=0; di<n && di<m-1; di++) {
    int k=0;
    for(int dj=di; dj>=0; dj--)
      out[di] += data[dj] * kernel[k++];
  }
 
  return(0);
}

corDispersion()

int corDispersion	(	double *	x,
		double *	y,
		const int	n,
		const double	tau,
		double *	tmp )

Correct time-activity curve for the effect of dispersion.

Data must be noise-free and have very short sampling intervals. The units of rate constants must be related to the TAC time units; 1/min and min, or 1/sec and sec.

Returns

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

Author

Vesa Oikonen

See also

simDispersion, simDelay

Parameters

x	Array of sample times; must be in increasing order and >=0.
y	Array of sample values, which will be replaced here by dispersion corrected values.
n	Nr of samples.
tau	Dispersion time constant (zero if no dispersion); in same time unit as sample times.
tmp	Array for temporary data, for at least n samples; enter NULL to let function to allocate and free the temporary space.

Definition at line 88 of file simdispersion.c.

  {
  /* Check input */
  if(x==NULL || y==NULL || n<2) return 1;
  if(tau<0.0) return 2;
  if(x[0]<0.0) return 3;
 
  /* Allocate memory if not allocated by user */
  double *buf;
  if(tmp!=NULL) buf=tmp; else buf=(double*)malloc(n*sizeof(double));
  if(buf==NULL) return 3;
 
  /* Integrate measured data */
  buf[0]=0.5*y[0]*x[0];
  for(int i=1; i<n; i++) buf[i]=buf[i-1]+0.5*(y[i]+y[i-1])*(x[i]-x[i-1]);
 
  /* Calculate the integral of dispersion corrected data */
  for(int i=0; i<n; i++) buf[i]+=tau*y[i];
 
  /* Calculate dispersion corrected curve from its integral */
  {
    int i=0;
    double dt=x[i];
    if(dt>0.0) y[i]=2.0*buf[i]/dt; else y[i]=0.0;
    for(i=1; i<n-1; i++) {
      dt=x[i+1]-x[i-1]; if(dt>0.0) y[i]=(buf[i+1]-buf[i-1])/dt; else y[i]=y[i-1];
    }
    dt=x[i]-x[i-1]; if(dt>0.0) y[i]=(buf[i]-buf[i-1])/dt; else y[i]=y[i-1];
  }
 
  if(tmp==NULL) free(buf);
  return 0;
}

icmparcAddMetabolites()

int icmparcAddMetabolites	(	ICMPARC *	d,
		unsigned const int	mNr )

Add sub-structures for metabolite(s) in parent ICMPARC structure.

See also

ICMPARC, icmparcInit, icmparcFree, icmparcAllocateTACs

Returns

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

Parameters

d	Pointer to parent ICMPARC; can itself be a sub-structure.
mNr	Number of metabolites to add.

Definition at line 43 of file simblood.c.

  {
  if(d==NULL) return(1);
  if(d->mNr>0 && d->metabolite!=NULL) { // delete any previous metabolite(s)
    for(unsigned int i=0; i<d->mNr; i++) icmparcFree(&d->metabolite[i]);
    free(d->metabolite); d->mNr=0; d->metabolite=(ICMPARC*)NULL;
  }
  if(mNr==0) return(0);
 
  d->metabolite=(ICMPARC*)malloc(sizeof(ICMPARC)*mNr);
  if(d->metabolite==NULL) return(3);
  d->mNr=mNr;
  for(unsigned int i=0; i<mNr; i++) {
    icmparcInit(&d->metabolite[i]);
    d->metabolite[i].parent=d;
  }
  return(0);
}

icmparcAllocateTACs()

int icmparcAllocateTACs	(	ICMPARC *	d,
		unsigned const int	sNr,
		const int	sub )

Allocate space for optional TACs inside ICMPARC structure.

See also

ICMPARC, icmparcInit, icmparcAddMetabolites, icmparcFree

Returns

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

Parameters

d	Pointer to ICMPARC.
sNr	Number of TAC samples (array lengths); set to zero to just free previous data.
sub	Allocate (1) or do not allocate (0) space for TACs inside sub-structures.

Definition at line 72 of file simblood.c.

  {
  if(d==NULL) return(1);
 
  /* Delete any existing TAC data */
  if(d->ic_BV!=NULL) free(d->ic_BV);
  if(d->ic_TS!=NULL) free(d->ic_TS);
  if(d->ic_TF!=NULL) free(d->ic_TF);
  d->ic_BV=d->ic_TS=d->ic_TF=(double*)NULL;
  if(d->c_BA!=NULL) free(d->c_BA);
  if(d->c_BV!=NULL) free(d->c_BV);
  if(d->c_TS!=NULL) free(d->c_TS);
  if(d->c_TF!=NULL) free(d->c_TF);
  d->c_BA=d->c_BV=d->c_TS=d->c_TF=(double*)NULL;
  if(sub!=0 && d->mNr>0 && d->metabolite!=NULL) {
    for(unsigned int i=0; i<d->mNr; i++) icmparcAllocateTACs(&d->metabolite[i], 0, sub);
  }
  if(sNr==0) return(0);
 
  /* Allocate space */
  d->c_BA=(double*)malloc(sNr*sizeof(double));
  d->c_BV=(double*)malloc(sNr*sizeof(double));
  d->c_TS=(double*)malloc(sNr*sizeof(double));
  d->c_TF=(double*)malloc(sNr*sizeof(double));
  d->ic_BV=(double*)malloc(sNr*sizeof(double));
  d->ic_TS=(double*)malloc(sNr*sizeof(double));
  d->ic_TF=(double*)malloc(sNr*sizeof(double));
  if(d->c_BA==NULL || d->c_BV==NULL || d->c_TS==NULL || d->c_TF==NULL ||
     d->ic_BV==NULL || d->ic_TS==NULL || d->ic_TF==NULL) {
    icmparcAllocateTACs(d, 0, sub); return(3);
  }
  if(sub!=0 && d->mNr>0 && d->metabolite!=NULL) {
    int ret=0;
    for(unsigned int i=0; i<d->mNr && !ret; i++)
      ret=icmparcAllocateTACs(&d->metabolite[i], sNr, sub);
    if(ret) return(ret);
  }
 
  return(0);
}

Referenced by icmparcAllocateTACs().

icmparcFree()

void icmparcFree

(

ICMPARC *

d

)

Free memory that is (optionally) allocated inside the ICMPARC structure or in its sub-structures. Field values are not changed, except for pointers to allocated memory that are set to NULL.

See also

ICMPARC, icmparcInit, icmparcAddMetabolites, icmparcAllocateTACs

Parameters

d	Pointer to ICMPARC.

Definition at line 125 of file simblood.c.

  {
  if(d==NULL) return;
  /* Free sub-structures recursively */
  if(d->mNr>0 && d->metabolite!=NULL) {
    for(unsigned int i=0; i<d->mNr; i++) icmparcFree(&d->metabolite[i]);
    free(d->metabolite); d->mNr=0; d->metabolite=(ICMPARC*)NULL;
  }
  /* Free optional TACs */
  if(d->ic_BV!=NULL) free(d->ic_BV);
  if(d->ic_TS!=NULL) free(d->ic_TS);
  if(d->ic_TF!=NULL) free(d->ic_TF);
  d->ic_BV=d->ic_TS=d->ic_TF=(double*)NULL;
  if(d->c_BA!=NULL) free(d->c_BA);
  if(d->c_BV!=NULL) free(d->c_BV);
  if(d->c_TS!=NULL) free(d->c_TS);
  if(d->c_TF!=NULL) free(d->c_TF);
  d->c_BA=d->c_BV=d->c_TS=d->c_TF=(double*)NULL;
}

Referenced by icmparcAddMetabolites(), and icmparcFree().

icmparcInit()

void icmparcInit

(

ICMPARC *

d

)

Initiate the ICMPARC structure before any use.

See also

ICMPARC, icmparcAddMetabolites, icmparcFree, icmparcAllocateTACs

Parameters

d	Pointer to ICMPARC.

Definition at line 20 of file simblood.c.

  {
  if(d==NULL) return;
  d->name[0]=(char)0;
  d->Ti=d->Tdur=d->Irate=0.0;
  d->k_BV_BA=0.0;
  d->k_BA_U=d->k_BA_TF=d->k_BA_TS=0.0;
  d->k_TF_BV=d->k_TS_BV=0.0;
  d->mNr=0; d->metabolite=(ICMPARC*)NULL;
  d->parent=(ICMPARC*)NULL;
  d->kp_BV=d->kp_TF=d->kp_TS=0.0;
  d->ic_BV=d->ic_TS=d->ic_TF=(double*)NULL;
  d->c_BA=d->c_BV=d->c_TS=d->c_TF=(double*)NULL;
}

Referenced by icmparcAddMetabolites().

simBM1()

int simBM1	(	double *	t,
		double *	ct,
		const int	nr,
		const double	k,
		double *	cp )

Calculate metabolite correct input curve from total input curve using simplistic 1K model dCmetabolite(t)/dt = k*Cparent(t) and since Ctotal(t) = Cparent(t) + Cmetabolite(t), we calculate metabolite corrected curve as dCparent(t)/dt = dCtotal(t)/dt - k*Cparent(t). Note that this model may result in negative parent values.

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

See also

simC1, simDispersion

Returns

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

Author

Vesa Oikonen

Parameters

t	Array of time values.
ct	Array of total activities.
nr	Number of values in TACs.
k	Rate constant of the model.
cp	Pointer for metabolite corrected data array to be simulated; must be allocated.

Definition at line 275 of file simblood.c.

  {
  int i;
  double dt2;
  double cp_last, t_last;
  double cpi, cpi_last;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(t==NULL || ct==NULL || cp==NULL) return 2;
 
  /* Check actual parameter number */
  if(!(k>=0.0)) return 3;
 
  /* Calculate curves */
  t_last=0.0; if(t[0]<t_last) t_last=t[0]; 
  cp_last=0.0; cpi=0.0; cpi_last=0.0; cpi_last=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) {
      cp[i] = (ct[i] - k*(cpi_last + dt2*cp_last)) / (1.0 + dt2*k);
      cpi = cpi_last + dt2*(cp_last+cp[i]);
    }
    /* prepare to the next loop */
    t_last=t[i]; cp_last=cp[i]; cpi_last=cpi;
  }
 
  return 0;
}

simBTAC()

int simBTAC	(	double *	t,
		const unsigned int	nr,
		ICMPARC *	p,
		double *	c_BA )

Simulate BTAC using compartmental model. NOT FUNCTIONAL!

See also

ICMPARC, icmparcInit, icmparcAddMetabolites, icmparcAllocateTACs, simDispersion, simSamples, parExampleTTACs, parExamplePerfectBolus, simDelay

Returns

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

Parameters

t	Array of sample time values; must be non-negative and in increasing order.
nr	Number of values (samples) in TACs.
p	Pointer to data structure containing simulation parameters.

See also

icmparcInit, icmparcAddMetabolites, icmparcAllocateTACs

Parameters

c_BA	Pointer for arterial BTAC array to be simulated; must be allocated. NULL can be given, if TAC space is allocated inside p structure.

Definition at line 154 of file simblood.c.

  {
  /* Check for data */
  if(nr<2 || t==NULL) return(1);
  if(p==NULL || (c_BA==NULL && p->c_BA==NULL)) return(2);
  if(!(t[0]>=0.0)) return(3);
  if(p->mNr>0 && p->metabolite==NULL) return(4);
 
  /* Calculate sum rates of metabolism */
  double km_BV=0.0, km_TS=0.0, km_TF=0.0;
  for(unsigned int i=0; i<p->mNr; i++) {
    if(p->metabolite[i].kp_BV>0.0) km_BV+=p->metabolite[i].kp_BV;
    if(p->metabolite[i].kp_TS>0.0) km_TS+=p->metabolite[i].kp_TS;
    if(p->metabolite[i].kp_TF>0.0) km_TF+=p->metabolite[i].kp_TF;
  }
 
  /* Calculate compartmental TACs */
  double t_last=0.0;
  double ic_BA_last=0.0, c_BA_last=0.0;
  double ic_BV_last=0.0, c_BV_last=0.0;
  double ic_TS_last=0.0, c_TS_last=0.0;
  double ic_TF_last=0.0, c_TF_last=0.0;
  for(unsigned int i=0; i<nr; i++) {
    double dt2=0.5*(t[i]-t_last); // delta time / 2
    if(!(dt2>=0.0)) return(3);
    else if(dt2==0.0) {
      if(c_BA!=NULL) c_BA[i]=c_BA_last;
      if(p->ic_BV!=NULL) p->ic_BV[i]=ic_BV_last;
      if(p->ic_TS!=NULL) p->ic_TS[i]=ic_TS_last;
      if(p->ic_TF!=NULL) p->ic_TF[i]=ic_TF_last;
      if(p->c_BA!=NULL) p->c_BA[i]=c_BA_last;
      if(p->c_BV!=NULL) p->c_BV[i]=c_BV_last;
      if(p->c_TS!=NULL) p->c_TS[i]=c_TS_last;
      if(p->c_TF!=NULL) p->c_TF[i]=c_TF_last;
      continue;
    }
    /* Infusion integral; before infusion start 0; during infusion: time*level; 
       thereafter (infusion time)*level and level is set to 1. */
    double ii=0.0;
    if(p->Ti>=0.0 && p->Tdur>0.0 && t[i]>p->Ti) {
      if(t[i] <= p->Ti + p->Tdur) ii=p->Irate*(t[i]-p->Ti); else ii=p->Irate*p->Tdur;
    }
    /* Helpers */
    double rBA=1.0/(1.0+dt2*(p->k_BA_TS+p->k_BA_TF+p->k_BA_U));
    double rTS=1.0/(1.0+dt2*(p->k_TS_BV+km_TS));
    double rTF=1.0/(1.0+dt2*(p->k_TF_BV+km_TF));
    double eBA=ic_BA_last+dt2*c_BA_last;
    double eBV=ic_BV_last+dt2*c_BV_last;
    double eTS=ic_TS_last+dt2*c_TS_last;
    double eTF=ic_TF_last+dt2*c_TF_last;
if(0) {
  printf("dt2=%g\n", dt2);
  printf("%g %g %g %g %g %g %g\n", rBA, rTS, rTF, eBA, eBV, eTS, eTF);
}
    /* Get integral of parent compound (if available) in compartments where it could be transformed 
       into current compound */
    double icp_BV=0.0, icp_TS=0.0, icp_TF=0.0;
    if(p->parent!=NULL) {
      // copy pre-calculated values
      if(p->parent->ic_BV!=NULL) icp_BV=p->parent->ic_BV[i];
      if(p->parent->ic_TS!=NULL) icp_TS=p->parent->ic_TS[i];
      if(p->parent->ic_TF!=NULL) icp_TF=p->parent->ic_TF[i];
    }
    /* Calculate integral of BV TAC */
    double ic_BV;
    ic_BV = eBV + dt2*ii - p->kp_BV*dt2*icp_BV 
            + dt2*dt2*rBA*eBA*(p->k_TS_BV*p->k_BA_TS*rTS + p->k_TF_BV*p->k_BA_TF*rTF)
            + p->k_TS_BV*dt2*rTS*(dt2*icp_TS + eTS) + p->k_TF_BV*dt2*rTF*(dt2*icp_TF + eTF);
    ic_BV /= 1.0 + dt2*(p->k_BV_BA + km_BV - dt2*dt2*rBA*p->k_BV_BA*(p->k_TS_BV*p->k_BA_TS*rTS + p->k_TF_BV*p->k_BA_TF*rTF));
    /* Calculate integral of BA TAC */
    double ic_BA = rBA * (dt2*p->k_BV_BA*ic_BV + eBA);
    /* Calculate integrals of TTACs */
    double ic_TS = rTS * (dt2*p->k_BA_TS*ic_BA + dt2*p->kp_TS*icp_TS + eTS);
    double ic_TF = rTF * (dt2*p->k_BA_TF*ic_BA + dt2*p->kp_TF*icp_TF + eTF);
    /* Prepare for the next sample */
    t_last=t[i];
    ic_BV_last=ic_BV; c_BV_last=(ic_BV-eBV)/dt2;
    ic_BA_last=ic_BA; c_BA_last=(ic_BA-eBA)/dt2;
    ic_TS_last=ic_TS; c_TS_last=(ic_TS-eTS)/dt2;
    ic_TF_last=ic_TF; c_TF_last=(ic_TF-eTF)/dt2;
    /* Copy sample concentration to output TAC */
    if(c_BA!=NULL) c_BA[i]=c_BA_last;
    /* If requested for metabolite calculations, copy also the compartmental TAC integrals */
    if(p->ic_BV!=NULL) p->ic_BV[i]=ic_BV;
    if(p->ic_TS!=NULL) p->ic_TS[i]=ic_TS;
    if(p->ic_TF!=NULL) p->ic_TF[i]=ic_TF;
    /* If requested, copy also the compartmental TACs */
    if(p->c_BA!=NULL) p->c_BA[i]=c_BA_last;
    if(p->c_BV!=NULL) p->c_BV[i]=c_BV_last;
    if(p->c_TS!=NULL) p->c_TS[i]=c_TS_last;
    if(p->c_TF!=NULL) p->c_TF[i]=c_TF_last;
  }
  return(0);
}

simC1()

int simC1	(	double *	t,
		double *	ca,
		const int	nr,
		const double	k1,
		const double	k2,
		double *	ct )

Simulate tissue TAC using 1 tissue compartmental model and plasma TAC, at plasma TAC times.

Memory for ct must be allocated in the calling program.

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

See also

simC1_d, simC1_i, simDispersion, simMBF, simC2, simOxygen, simC3vs, convolve1D, simC1DI, simSamples, parExampleTTACs, parExamplePerfectBolus

Returns

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

Author

Vesa Oikonen

Parameters

t	Array of time values.
ca	Array of arterial activities.
nr	Number of values in TACs.
k1	Rate constant of the model.
k2	Rate constant of the model.
ct	Pointer for TAC array to be simulated; must be allocated.

Definition at line 93 of file sim1cm.c.

  {
  int i;
  double dt2;
  double cai, ca_last, t_last;
  double ct1, ct1_last;
  double ct1i, ct1i_last;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(t==NULL || ca==NULL || ct==NULL) return 2;
 
  /* Check actual parameter number */
  if(!(k1>=0.0)) return 3;
 
  /* Calculate curves */
  t_last=0.0; if(t[0]<t_last) t_last=t[0]; 
  cai=ca_last=0.0;
  ct1_last=ct1i_last=0.0;
  ct1=ct1i=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 integral */
      cai+=(ca[i]+ca_last)*dt2;
      /* tissue compartment and its integral */
      ct1 = (k1*cai - k2*(ct1i_last+dt2*ct1_last)) / (1.0 + dt2*k2);
      ct1i = ct1i_last + dt2*(ct1_last+ct1);
    }
    /* copy values to argument arrays */
    ct[i]=ct1;
    /* prepare to the next loop */
    t_last=t[i]; ca_last=ca[i];
    ct1_last=ct1; ct1i_last=ct1i;
  }
 
  return 0;
}

Referenced by bfm1TCM(), and simDispersion().

simC1_d()

int simC1_d	(	double *	t,
		double *	ca,
		const int	nr,
		const double	k1,
		const double	k2,
		double *	ct )

Simulate tissue TAC using 1 tissue compartmental model and input TAC, at input TAC times.

This version is directly based on ODE, and does not use integral of input TAC.

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

See also

simC1, simDelay, simDispersion

Returns

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

Author

Vesa Oikonen

Parameters

t	Array of time values.
ca	Array of arterial activities.
nr	Number of values in TACs.
k1	Rate constant of the model.
k2	Rate constant of the model.
ct	Pointer for TAC array to be simulated; must be allocated.

Definition at line 229 of file sim1cm.c.

  {
  int i;
  double dt2;
  double cap, ctp, tp; // 'p' marks previous
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(t==NULL || ca==NULL || ct==NULL) return 2;
 
  /* Check actual parameter number */
  if(!(k1>=0.0)) return 3;
 
  /* Set 'previous' sample */
  tp=0.0; if(t[0]<tp) tp=t[0]; 
  cap=ctp=0.0;
  /* Calculate curves */
  for(i=0; i<nr; i++) {
    /* delta time / 2 */
    dt2=0.5*(t[i]-tp);
    /* calculate values */
    if(dt2<0.0) {
      return 5;
    } else if(dt2>=0.0) { // must be >= because ct[0] is not set elsewhere
      ct[i]=(k1*dt2*(cap+ca[i]) + (1.0 - k2*dt2)*ctp) / (1.0 + k2*dt2);
    }
    /* prepare to the next sample */
    tp=t[i]; cap=ca[i];
    ctp=ct[i];
  }
 
  return 0;
}

simC1_i()

int simC1_i	(	double *	t,
		double *	cai,
		const int	nr,
		const double	k1,
		const double	k2,
		double *	ct )

Simulate tissue TAC using 1 tissue compartmental model and plasma TAC, at plasma TAC times.

This version uses integral of arterial TAC as input function. Only advantage over simC1() is that the calculation of integral can be fully controlled and possibly more precise in some situations.

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

See also

simC1, simC2_i, simDelay, simDispersion

Returns

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

Author

Vesa Oikonen

Parameters

t	Array of time values.
cai	Array of AUC 0-t of arterial activities.
nr	Number of values in TACs.
k1	Rate constant of the model.
k2	Rate constant of the model.
ct	Pointer for TAC array to be simulated; must be allocated.

Definition at line 164 of file sim1cm.c.

  {
  int i;
  double dt2;
  double t_last;
  double ct1, ct1_last;
  double ct1i, ct1i_last;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(t==NULL || cai==NULL || ct==NULL) return 2;
 
  /* Check actual parameter number */
  if(!(k1>=0.0)) return 3;
 
  /* Calculate curves */
  t_last=0.0; if(t[0]<t_last) t_last=t[0]; 
  ct1_last=ct1i_last=0.0;
  ct1=ct1i=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) {
      /* tissue compartment and its integral */
      ct1 = (k1*cai[i] - k2*(ct1i_last+dt2*ct1_last)) / (1.0 + dt2*k2);
      ct1i = ct1i_last + dt2*(ct1_last+ct1);
    }
    /* copy values to argument arrays */
    ct[i]=ct1;
    /* prepare to the next loop */
    t_last=t[i]; ct1_last=ct1; ct1i_last=ct1i;
  }
 
  return 0;
}

Referenced by bfm1TCM(), and bfmSRTM().

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);
}

simC2()

int simC2	(	double *	t,
		double *	ca,
		const int	nr,
		const double	k1,
		const double	k2,
		const double	k3,
		const double	k4,
		double *	ct,
		double *	cta,
		double *	ctb )

Simulate tissue TAC using two-tissue compartment model and plasma TAC, at plasma TAC times.

Memory for ct must be allocated in the calling program. To retrieve the separate tissue compartment TACs, pointer to allocated memory for cta and/or ctb can be given; if compartmental TACs are not required, NULL pointer can be given instead.

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

simC2_i, simC1, simC1DI, simC2l, simC3s, simC3p, simSRTM, simSamples, simDispersion, parExampleTTACs, parExamplePerfectBolus

Parameters

t	Array of time values.
ca	Array of arterial activities.
nr	Number of values in TACs.
k1	Rate constant of the model.
k2	Rate constant of the model.
k3	Rate constant of the model.
k4	Rate constant of the model.
ct	Pointer for TAC array to be simulated; must be allocated.
cta	Pointer for 1st compartment TAC to be simulated, or NULL.
ctb	Pointer for 2nd compartment TAC to be simulated, or NULL.

Definition at line 31 of file sim2cm.c.

  {
  int i;
  double dt2, r, u, v;
  double cai, ca_last, t_last;
  double ct1, ct1_last, ct2, ct2_last;
  double ct1i, ct1i_last, ct2i, ct2i_last;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(t==NULL || ca==NULL || ct==NULL) return 2;
 
  /* Check parameters */
  if(k1<0.0) return 3;
 
  /* Calculate curves */
  t_last=0.0; if(t[0]<t_last) t_last=t[0];
  cai=ca_last=0.0;
  ct1_last=ct2_last=ct1i_last=ct2i_last=0.0;
  ct1=ct2=ct1i=ct2i=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 integral */
      cai+=(ca[i]+ca_last)*dt2;
      /* Calculate partial results */
      r=1.0+k4*dt2;
      u=ct1i_last+dt2*ct1_last;
      v=ct2i_last+dt2*ct2_last;
      /* 1st tissue compartment and its integral */
      ct1 = ( k1*cai - (k2 + (k3/r))*u + (k4/r)*v ) / ( 1.0 + dt2*(k2 + (k3/r)) );
      ct1i = ct1i_last + dt2*(ct1_last+ct1);
      /* 2nd tissue compartment and its integral */
      ct2 = (k3*ct1i - k4*v) / r;
      ct2i = ct2i_last + dt2*(ct2_last+ct2);
    }
    /* copy values to argument arrays */
    ct[i]=ct1+ct2;
    if(cta!=NULL) cta[i]=ct1;
    if(ctb!=NULL) ctb[i]=ct2;
    /* prepare to the next loop */
    t_last=t[i]; ca_last=ca[i];
    ct1_last=ct1; ct1i_last=ct1i;
    ct2_last=ct2; ct2i_last=ct2i;
  }
 
  return 0;
}

simC2_i()

int simC2_i	(	double *	t,
		double *	cai,
		const int	nr,
		const double	k1,
		const double	k2,
		const double	k3,
		const double	k4,
		double *	ct,
		double *	cta,
		double *	ctb )

Simulate tissue TAC using two-tissue compartment model and plasma TAC, at plasma TAC times.

Memory for ct must be allocated in the calling program. To retrieve the separate tissue compartment TACs, pointer to allocated memory for cta and/or ctb can be given; if compartmental TACs are not required, NULL pointer can be given instead.

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

This version uses integral of arterial TAC as input function. Only advantage over simC2() is that the calculation of integral can be fully controlled and possibly more precise in some situations.

Returns

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

Author

Vesa Oikonen

See also

simC2, simC1, simC1_i, simC3s, simC3p, simSamples, simDelay

Parameters

t	Array of time values.
cai	Array of AUC 0-t of arterial activities.
nr	Number of values in TACs.
k1	Rate constant of the model.
k2	Rate constant of the model.
k3	Rate constant of the model.
k4	Rate constant of the model.
ct	Pointer for TAC array to be simulated; must be allocated.
cta	Pointer for 1st compartment TAC to be simulated, or NULL.
ctb	Pointer for 2nd compartment TAC to be simulated, or NULL.

Definition at line 124 of file sim2cm.c.

  {
  int i;
  double dt2, r, u, v;
  double t_last;
  double ct1, ct1_last, ct2, ct2_last;
  double ct1i, ct1i_last, ct2i, ct2i_last;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(t==NULL || cai==NULL || ct==NULL) return 2;
 
  /* Check parameters */
  if(k1<0.0) return 3;
 
  /* Calculate curves */
  t_last=0.0; if(t[0]<t_last) t_last=t[0];
  ct1_last=ct2_last=ct1i_last=ct2i_last=0.0;
  ct1=ct2=ct1i=ct2i=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) {
      /* Calculate partial results */
      r=1.0+k4*dt2;
      u=ct1i_last+dt2*ct1_last;
      v=ct2i_last+dt2*ct2_last;
      /* 1st tissue compartment and its integral */
      ct1 = ( k1*cai[i] - (k2 + (k3/r))*u + (k4/r)*v ) / ( 1.0 + dt2*(k2 + (k3/r)) );
      ct1i = ct1i_last + dt2*(ct1_last+ct1);
      /* 2nd tissue compartment and its integral */
      ct2 = (k3*ct1i - k4*v) / r;
      ct2i = ct2i_last + dt2*(ct2_last+ct2);
    }
    /* copy values to argument arrays */
    ct[i]=ct1+ct2;
    if(cta!=NULL) cta[i]=ct1;
    if(ctb!=NULL) ctb[i]=ct2;
    /* prepare to the next loop */
    t_last=t[i];
    ct1_last=ct1; ct1i_last=ct1i;
    ct2_last=ct2; ct2i_last=ct2i;
  }
 
  return 0;
}

simC2l()

int simC2l	(	double *	t,
		double *	ca,
		const int	nr,
		const double	k1,
		const double	k2,
		const double	k3,
		const double	kLoss,
		double *	ct,
		double *	cta,
		double *	ctb )

Simulate tissue TAC using 2 tissue compartment model (in series) and plasma TAC, at plasma TAC times. In contrary to the common model, kLoss represents a direct loss rate from the 2nd tissue compartment to venous plasma.

Memory for ct must be allocated in the calling program. To retrieve the separate tissue compartment TACs, pointer to allocated memory for cta and ctb can be given; if compartmental TACs are not required, NULL pointer can be given instead.

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

See also

simC2, simC2vl, simC3vs, simC3vpKLoss

Returns

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

Author

Vesa Oikonen

Parameters

t	Array of time values.
ca	Array of arterial activities.
nr	Number of values in TACs.
k1	Rate constant of the model.
k2	Rate constant of the model.
k3	Rate constant of the model.
kLoss	Rate constant of the model.
ct	Pointer for TAC array to be simulated; must be allocated.
cta	Pointer for 1st compartment TAC to be simulated, or NULL.
ctb	Pointer for 2nd compartment TAC to be simulated, or NULL.

Definition at line 169 of file simkloss.c.

  {
  int i;
  double b, c, dt2;
  double cai, ca_last, t_last;
  double ct1, ct1_last, ct2, ct2_last;
  double ct1i, ct1i_last, ct2i, ct2i_last;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(ct==NULL) return 2;
 
  /* Check actual parameter number */
  if(k1<0.0) return 3;
 
  /* Calculate curves */
  t_last=0.0; if(t[0]<t_last) t_last=t[0];
  cai=ca_last=0.0;
  ct1_last=ct2_last=ct1i_last=ct2i_last=ct1=ct2=ct1i=ct2i=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 integral */
      cai+=(ca[i]+ca_last)*dt2;
      /* partial results */
      b=ct1i_last+dt2*ct1_last;
      c=ct2i_last+dt2*ct2_last;
      /* 1st tissue compartment and its integral */
      ct1 = (k1*cai - (k2+k3)*b) / (1.0 + (k2+k3)*dt2 );
      ct1i = ct1i_last + dt2*(ct1_last+ct1);
      /* 2nd tissue compartment and its integral */
      ct2 = (k3*ct1i - kLoss*c) / (1.0 + kLoss*dt2);
      ct2i = ct2i_last + dt2*(ct2_last+ct2);
    }
    /* copy values to argument arrays */
    ct[i]=ct1+ct2;
    if(cta!=NULL) cta[i]=ct1;
    if(ctb!=NULL) ctb[i]=ct2;
    /* prepare to the next loop */
    t_last=t[i]; ca_last=ca[i];
    ct1_last=ct1; ct1i_last=ct1i;
    ct2_last=ct2; ct2i_last=ct2i;
  }
 
  return 0;
}

simC2vl()

int simC2vl	(	double *	t,
		double *	ca,
		double *	cb,
		const int	nr,
		const double	k1,
		const double	k2,
		const double	k3,
		const double	kL,
		const double	f,
		const double	vb,
		const double	fa,
		const int	vvm,
		double *	cpet,
		double *	cta,
		double *	ctb,
		double *	ctab,
		double *	ctvb )

Simulate tissue TAC using 2 tissue compartment model and plasma TAC, at plasma TAC times, considering also arterial and venous vasculature. The efflux from 2nd tissue compartment (at rate kL) goes directly to blood.

Memory for cpet must be allocated in the calling program. To retrieve the separate tissue compartment TACs, pointer to allocated memory for cta, ctb, ctab, and/or ctvb can be given; if compartmental TACs are not required, NULL pointer can be given instead.

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.

See also

simC2, simC2l, simC3vs, simC3vpKLoss, simC3DIvs

Returns

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

Author

Vesa Oikonen

Parameters

t	Array of time values.
ca	Array of arterial plasma activities.
cb	Array of arterial blood activities.
nr	Number of values in TACs.
k1	Rate constant of the model.
k2	Rate constant of the model.
k3	Rate constant of the model.
kL	Rate constant of the model.
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.
cpet	Pointer for TAC array to be simulated; must be allocated.
cta	Pointer for 1st compartment TAC to be simulated, or NULL.
ctb	Pointer for 2nd compartment TAC to be simulated, or NULL.
ctab	Pointer for arterial TAC in tissue, or NULL.
ctvb	Pointer for venous TAC in tissue, or NULL.

Definition at line 261 of file simkloss.c.

  {
  int i;
  double dt2, b, c, va, vv;
  double cai, ca_last, t_last, dct, cvb;
  double ct1, ct1_last, ct2, ct2_last;
  double ct1i, ct1i_last, ct2i, ct2i_last;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(cpet==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];
  cai=ca_last=0.0;
  ct1_last=ct2_last=ct1i_last=ct2i_last=ct1=ct2=ct1i=ct2i=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 integral */
      cai+=(ca[i]+ca_last)*dt2;
      /* Calculate partial results */
      b=ct1i_last+dt2*ct1_last;
      c=ct2i_last+dt2*ct2_last;
      /* 1st tissue compartment and its integral */
      ct1 = (k1*cai - (k2+k3)*b) / (1.0 + (k2+k3)*dt2 );
      ct1i = ct1i_last + dt2*(ct1_last+ct1);
      /* 2nd tissue compartment and its integral */
      ct2 = (k3*ct1i - kL*c) / (1.0 + kL*dt2);
      ct2i = ct2i_last + dt2*(ct2_last+ct2);
    }
    /* Venous curve */
    if(f>0.) {dct = k1*ca[i] - k2*ct1 - kL*ct2; cvb = cb[i] - dct/f;}
    else cvb=cb[i];
 
    /* copy values to argument arrays */
    if(vvm==0) cpet[i]= va*cb[i] + vv*cvb + (1.0-vb)*(ct1+ct2);
    else cpet[i]= va*cb[i] + vv*cvb + (ct1+ct2);
    if(cta!=NULL) {cta[i]=ct1; if(vvm==0) cta[i]*=(1.0-vb);}
    if(ctb!=NULL) {ctb[i]=ct2; if(vvm==0) ctb[i]*=(1.0-vb);}
    if(ctab!=NULL) {ctab[i]=va*cb[i];}
    if(ctvb!=NULL) {ctvb[i]=vv*cvb;}
    /* prepare to the next loop */
    t_last=t[i]; ca_last=ca[i];
    ct1_last=ct1; ct1i_last=ct1i;
    ct2_last=ct2; ct2i_last=ct2i;
  }
 
  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;
}

simC3p()

int simC3p	(	double *	t,
		double *	ca,
		const int	nr,
		const double	k1,
		const double	k2,
		const double	k3,
		const double	k4,
		const double	k5,
		const double	k6,
		double *	ct,
		double *	cta,
		double *	ctb,
		double *	ctc )

Simulate tissue TAC using 1-3 tissue compartment model (2nd and 3rd compartments in parallel) and plasma TAC, at plasma TAC times.

Memory for ct must be allocated in the calling program. To retrieve the separate tissue compartment TACs, pointer to allocated memory for cta, ctb and/or ctc can be given; if compartmental TACs are not required, NULL pointer can be given instead.

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

simC3s, simC1, simC1DI, simC2, simC2l, simC3vp, simC3vs, simSamples, parExampleTTACs, parExamplePerfectBolus

Parameters

t	Array of time values.
ca	Array of arterial activities.
nr	Number of values in TACs.
k1	Rate constant of the model.
k2	Rate constant of the model.
k3	Rate constant of the model.
k4	Rate constant of the model.
k5	Rate constant of the model.
k6	Rate constant of the model.
ct	Pointer for TAC array to be simulated; must be allocated.
cta	Pointer for 1st compartment TAC to be simulated, or NULL.
ctb	Pointer for 2nd compartment TAC to be simulated, or NULL.
ctc	Pointer for 3rd compartment TAC to be simulated, or NULL.

Definition at line 33 of file sim3cmp.c.

  {
  int i;
  double dt2, r, s, u, v, w;
  double cai, ca_last, t_last;
  double ct1, ct1_last, ct2, ct2_last, ct3, ct3_last;
  double ct1i, ct1i_last, ct2i, ct2i_last, ct3i, ct3i_last;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(ct==NULL) return 2;
 
  /* Check parameters */
  if(k1<0.0) return 3;
 
  /* Calculate curves */
  t_last=0.0; if(t[0]<t_last) t_last=t[0];
  cai=ca_last=0.0;
  ct1_last=ct2_last=ct3_last=ct1i_last=ct2i_last=ct3i_last=0.0;
  ct1=ct2=ct3=ct1i=ct2i=ct3i=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 integral */
      cai+=(ca[i]+ca_last)*dt2;
      /* Calculate partial results */
      r=1.0+k4*dt2;
      s=1.0+k6*dt2;
      u=ct1i_last+dt2*ct1_last;
      v=ct2i_last+dt2*ct2_last;
      w=ct3i_last+dt2*ct3_last;
      /* 1st tissue compartment and its integral */
      ct1 = ( k1*cai - (k2 + (k3/r) + (k5/s))*u + (k4/r)*v + (k6/s)*w )
            / ( 1.0 + dt2*(k2 + (k3/r) + (k5/s)) );
      ct1i = ct1i_last + dt2*(ct1_last+ct1);
      /* 2nd tissue compartment and its integral */
      ct2 = (k3*ct1i - k4*v) / r;
      ct2i = ct2i_last + dt2*(ct2_last+ct2);
      /* 3rd tissue compartment and its integral */
      ct3 = (k5*ct1i - k6*w) / s;
      ct3i = ct3i_last + dt2*(ct3_last+ct3);
    }
    /* copy values to argument arrays; set very small values to zero */
    ct[i]=ct1+ct2+ct3;
    if(cta!=NULL) cta[i]=ct1;
    if(ctb!=NULL) ctb[i]=ct2;
    if(ctc!=NULL) ctc[i]=ct3;
    /* prepare to the next loop */
    t_last=t[i]; ca_last=ca[i];
    ct1_last=ct1; ct1i_last=ct1i;
    ct2_last=ct2; ct2i_last=ct2i;
    ct3_last=ct3; ct3i_last=ct3i;
  }
 
  return 0;
}

simC3s()

int simC3s	(	double *	t,
		double *	ca,
		const int	nr,
		double	k1,
		double	k2,
		double	k3,
		double	k4,
		double	k5,
		double	k6,
		double *	ct,
		double *	cta,
		double *	ctb,
		double *	ctc )

Simulate tissue TAC using 1-3 tissue compartment model (compartments in series) and plasma TAC, at plasma TAC times.

Memory for ct must be allocated in the calling program. To retrieve the separate tissue compartment TACs, pointer to allocated memory for cta, ctb and/or ctc can be given; if compartmental TACs are not required, NULL pointer can be given instead.

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

simC3vs, simC1, simC3p, simC3vp, simC2l, simDispersion, simDelay

Parameters

t	Array of time values.
ca	Array of arterial activities.
nr	Number of values in TACs.
k1	Rate constant of the model.
k2	Rate constant of the model.
k3	Rate constant of the model.
k4	Rate constant of the model.
k5	Rate constant of the model.
k6	Rate constant of the model.
ct	Pointer for TAC array to be simulated; must be allocated.
cta	Pointer for 1st compartment TAC to be simulated, or NULL.
ctb	Pointer for 2nd compartment TAC to be simulated, or NULL.
ctc	Pointer for 3rd compartment TAC to be simulated, or NULL.

Definition at line 32 of file sim3cms.c.

  {
  int i;
  double b, c, d, w, z, dt2;
  double cai, ca_last, t_last;
  double ct1, ct1_last, ct2, ct2_last, ct3, ct3_last;
  double ct1i, ct1i_last, ct2i, ct2i_last, ct3i, ct3i_last;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(ct==NULL) return 2;
 
  /* Check actual parameter number */
  if(k1<0.0) return 3;
  if(k3<=0.0) {k3=0.0; if(k2<=0.0) k2=0.0;}
  else if(k5<=0.0) {k5=0.0; if(k4<=0.0) k4=0.0;}
  else {if(k6<=0.0) k6=0.0;}
 
  /* Calculate curves */
  t_last=0.0; if(t[0]<t_last) t_last=t[0]; 
  cai=ca_last=0.0;
  ct1_last=ct2_last=ct3_last=ct1i_last=ct2i_last=ct3i_last=0.0;
  ct1=ct2=ct3=ct1i=ct2i=ct3i=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 integral */
      cai+=(ca[i]+ca_last)*dt2;
      /* partial results */
      b=ct1i_last+dt2*ct1_last;
      c=ct2i_last+dt2*ct2_last;
      d=ct3i_last+dt2*ct3_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*cai + (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);
      /* 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);
    }
    /* copy values to argument arrays */
    ct[i]=ct1+ct2+ct3;
    if(cta!=NULL) cta[i]=ct1;
    if(ctb!=NULL) ctb[i]=ct2;
    if(ctc!=NULL) ctc[i]=ct3;
    /* prepare to the next loop */
    t_last=t[i]; ca_last=ca[i];
    ct1_last=ct1; ct1i_last=ct1i;
    ct2_last=ct2; ct2i_last=ct2i;
    ct3_last=ct3; ct3i_last=ct3i;
  }
 
  return 0;
}

simC3vp()

int simC3vp	(	double *	t,
		double *	ca,
		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	f,
		const double	vb,
		const double	fa,
		const int	vvm,
		double *	cpet,
		double *	cta,
		double *	ctb,
		double *	ctc,
		double *	ctab,
		double *	ctvb )

Simulate tissue TAC using 1-3 tissue compartment model (2nd and 3rd compartments in parallel) and plasma TAC, at plasma TAC times, considering also arterial and venous vasculature.

Memory for ct must be allocated in the calling program. To retrieve the separate tissue compartment TACs, pointer to allocated memory for cta, ctb, ctc, ctab, and/or ctvb can be given; if compartmental TACs are not required, NULL pointer can be given instead.

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

simC3vs, simC3p, simC3s, simC3vpKLoss, simC4DIvp, simSamples, simDispersion, simDelay, parExampleTTACs, parExamplePerfectBolus

Parameters

t	Array of time values.
ca	Array of arterial plasma activities.
cb	Array of arterial blood activities.
nr	Number of values in TACs.
k1	Rate constant of the model.
k2	Rate constant of the model.
k3	Rate constant of the model.
k4	Rate constant of the model.
k5	Rate constant of the model.
k6	Rate constant of the model.
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.
cpet	Pointer for TAC array to be simulated; must be allocated.
cta	Pointer for 1st compartment TAC to be simulated, or NULL.
ctb	Pointer for 2nd compartment TAC to be simulated, or NULL.
ctc	Pointer for 3rd compartment TAC to be simulated, or NULL.
ctab	Pointer for arterial TAC in tissue, or NULL.
ctvb	Pointer for venous TAC in tissue, or NULL.

Definition at line 140 of file sim3cmp.c.

  {
  int i;
  double dt2, r, s, u, v, w, va, vv;
  double cai, ca_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;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(cpet==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];
  cai=ca_last=0.0;
  ct1_last=ct2_last=ct3_last=ct1i_last=ct2i_last=ct3i_last=0.0;
  ct1=ct2=ct3=ct1i=ct2i=ct3i=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 integral */
      cai+=(ca[i]+ca_last)*dt2;
      /* Calculate partial results */
      r=1.0+k4*dt2;
      s=1.0+k6*dt2;
      u=ct1i_last+dt2*ct1_last;
      v=ct2i_last+dt2*ct2_last;
      w=ct3i_last+dt2*ct3_last;
      /* 1st tissue compartment and its integral */
      ct1 = ( k1*cai - (k2 + (k3/r) + (k5/s))*u + (k4/r)*v + (k6/s)*w )
            / ( 1.0 + dt2*(k2 + (k3/r) + (k5/s)) );
      ct1i = ct1i_last + dt2*(ct1_last+ct1);
      /* 2nd tissue compartment and its integral */
      ct2 = (k3*ct1i - k4*v) / r;
      ct2i = ct2i_last + dt2*(ct2_last+ct2);
      /* 3rd tissue compartment and its integral */
      ct3 = (k5*ct1i - k6*w) / s;
      ct3i = ct3i_last + dt2*(ct3_last+ct3);
    }
    /* Venous curve */
    if(f>0.) {dct = k1*ca[i] - k2*ct1; cvb = cb[i] - dct/f;} else cvb=cb[i];
    /* copy values to argument arrays */
    if(vvm==0) cpet[i]= va*cb[i] + vv*cvb + (1.0-vb)*(ct1+ct2+ct3);
    else cpet[i]= va*cb[i] + vv*cvb + (ct1+ct2+ct3);
    if(cta!=NULL) {cta[i]=ct1; if(vvm==0) cta[i]*=(1.0-vb);}
    if(ctb!=NULL) {ctb[i]=ct2; if(vvm==0) ctb[i]*=(1.0-vb);}
    if(ctc!=NULL) {ctc[i]=ct3; if(vvm==0) ctc[i]*=(1.0-vb);}
    if(ctab!=NULL) {ctab[i]=va*cb[i];}
    if(ctvb!=NULL) {ctvb[i]=vv*cvb;}
    /* prepare to the next loop */
    t_last=t[i]; ca_last=ca[i];
    ct1_last=ct1; ct1i_last=ct1i;
    ct2_last=ct2; ct2i_last=ct2i;
    ct3_last=ct3; ct3i_last=ct3i;
  }
 
  return 0;
}

simC3vpKLoss()

int simC3vpKLoss	(	double *	t,
		double *	ca,
		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	kLoss,
		const double	f,
		const double	vb,
		const double	fa,
		const int	vvm,
		double *	cpet,
		double *	cta,
		double *	ctb,
		double *	ctc,
		double *	ctab,
		double *	ctvb )

Simulate tissue TAC using 3 tissue compartmental model with two parallel compartments, and plasma TAC, at plasma TAC sample times, considering also arterial and venous vasculature. The efflux from 3rd tissue compartment (C) goes directly to blood at rate kLoss.

Memory for cpet must be allocated in the calling program. To retrieve the separate tissue compartment TACs, pointer to allocated memory for cta, ctb, ctab, and/or ctvb can be given; if compartmental TACs are not required, NULL pointer can be given instead.

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.,

See also

simC2vl, simC2l, simC3vs, simC4DIvp, parExamplePerfectBolus

Returns

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

Author

Vesa Oikonen

Parameters

t	Array of sample times.
ca	Array of arterial plasma activities.
cb	Array of arterial blood activities.
nr	Number of sample values in TACs.
k1	Rate constant of the model.
k2	Rate constant of the model.
k3	Rate constant of the model.
k4	Rate constant of the model.
k5	Rate constant of the model.
k6	Rate constant of the model.
kLoss	Rate constant of the model.
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.
cpet	Pointer for TAC array to be simulated; must be allocated.
cta	Pointer for 1st tissue compartment TAC to be simulated, or NULL.
ctb	Pointer for 2nd compartment TAC to be simulated, or NULL.
ctc	Pointer for 3rd compartment TAC to be simulated, or NULL.
ctab	Pointer for arterial TAC in tissue, or NULL.
ctvb	Pointer for venous TAC in tissue, or NULL.

Definition at line 36 of file simkloss.c.

  {
  int i;
  double dt2, b, c, d, w, z, u, va, vv;
  double cai, ca_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;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(cpet==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];
  cai=ca_last=0.0;
  ct1_last=ct2_last=ct3_last=ct1i_last=ct2i_last=ct3i_last=0.0;
  ct1=ct2=ct3=ct1i=ct2i=ct3i=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 integral */
      cai+=(ca[i]+ca_last)*dt2;
      /* Calculate partial results */
      w = 1.0 + k4*dt2;
      z = 1.0 + dt2*(k6 + kLoss); 
      u = k2 + k3 + k5 - k3*k4*dt2/w - k5*k6*dt2/z;
      b = ct1i_last+dt2*ct1_last;
      c = ct2i_last+dt2*ct2_last;
      d = ct3i_last+dt2*ct3_last;
      /* 1st tissue compartment and its integral */
      ct1 = ( k1*cai - u*b + k4*c/w + k6*d/z ) / ( 1.0 + dt2*u);
      ct1i = ct1i_last + dt2*(ct1_last+ct1);
      /* 2nd tissue compartment and its integral */
      ct2 = (k3*ct1i - k4*c) / w;
      ct2i = ct2i_last + dt2*(ct2_last+ct2);
      /* 3rd tissue compartment and its integral */
      ct3 = (k5*ct1i - (k6 + kLoss)*d) / z;
      ct3i = ct3i_last + dt2*(ct3_last+ct3);
    }
    /* Venous curve */
    if(f>0.) {dct = k1*ca[i] - k2*ct1 - kLoss*ct3; cvb = cb[i] - dct/f;}
    else cvb=cb[i];
    /* copy values to argument arrays */
    if(vvm==0) cpet[i]= va*cb[i] + vv*cvb + (1.0-vb)*(ct1+ct2+ct3);
    else cpet[i]= va*cb[i] + vv*cvb + (ct1+ct2+ct3);
    if(cta!=NULL) {cta[i]=ct1; if(vvm==0) cta[i]*=(1.0-vb);}
    if(ctb!=NULL) {ctb[i]=ct2; if(vvm==0) ctb[i]*=(1.0-vb);}
    if(ctc!=NULL) {ctc[i]=ct3; if(vvm==0) ctc[i]*=(1.0-vb);}
    if(ctab!=NULL) {ctab[i]=va*cb[i];}
    if(ctvb!=NULL) {ctvb[i]=vv*cvb;}
    /* prepare to the next loop */
    t_last=t[i]; ca_last=ca[i];
    ct1_last=ct1; ct1i_last=ct1i;
    ct2_last=ct2; ct2i_last=ct2i;
    ct3_last=ct3; ct3i_last=ct3i;
  }
  
  return 0;
}

simC3vs()

int simC3vs	(	double *	t,
		double *	ca,
		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	f,
		const double	vb,
		const double	fa,
		const int	vvm,
		double *	cpet,
		double *	cta,
		double *	ctb,
		double *	ctc,
		double *	ctab,
		double *	ctvb )

Simulate tissue TAC using 1-3 tissue compartment model (compartments in series) and plasma TAC, at plasma TAC times, considering also arterial and venous vasculature.

Memory for ct must be allocated in the calling program. To retrieve the separate tissue compartment TACs, pointer to allocated memory for cta, ctb, ctc, ctab, and/or ctvb can be given; if compartmental TACs are not required, NULL pointer can be given instead.

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

simC3vp, simC3s, simC3p, simC2l, simC3vpKLoss, simC3DIvs, simC4DIvs, parExamplePerfectBolus

Parameters

t	Array of time values.
ca	Array of arterial plasma activities.
cb	Array of arterial blood activities.
nr	Number of values in TACs.
k1	Rate constant of the model.
k2	Rate constant of the model.
k3	Rate constant of the model.
k4	Rate constant of the model.
k5	Rate constant of the model.
k6	Rate constant of the model.
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.
cpet	Pointer for TAC array to be simulated; must be allocated.
cta	Pointer for 1st compartment TAC to be simulated, or NULL.
ctb	Pointer for 2nd compartment TAC to be simulated, or NULL.
ctc	Pointer for 3rd compartment TAC to be simulated, or NULL.
ctab	Pointer for arterial TAC in tissue, or NULL.
ctvb	Pointer for venous TAC in tissue, or NULL.

Definition at line 143 of file sim3cms.c.

  {
  int i;
  double b, c, d, w, z, dt2, va, vv;
  double cai, ca_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;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(cpet==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];
  cai=ca_last=0.0;
  ct1_last=ct2_last=ct3_last=ct1i_last=ct2i_last=ct3i_last=0.0;
  ct1=ct2=ct3=ct1i=ct2i=ct3i=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 integral */
      cai+=(ca[i]+ca_last)*dt2;
      /* partial results */
      b=ct1i_last+dt2*ct1_last;
      c=ct2i_last+dt2*ct2_last;
      d=ct3i_last+dt2*ct3_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*cai + (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);
      /* 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*ca[i] - k2*ct1; cvb = cb[i] - dct/f;} else cvb=cb[i];
    /* copy values to argument arrays */
    if(vvm==0) cpet[i]= va*cb[i] + vv*cvb + (1.0-vb)*(ct1+ct2+ct3);
    else cpet[i]= va*cb[i] + vv*cvb + (ct1+ct2+ct3);
    if(cta!=NULL) {cta[i]=ct1; if(vvm==0) cta[i]*=(1.0-vb);}
    if(ctb!=NULL) {ctb[i]=ct2; if(vvm==0) ctb[i]*=(1.0-vb);}
    if(ctc!=NULL) {ctc[i]=ct3; if(vvm==0) ctc[i]*=(1.0-vb);}
    if(ctab!=NULL) {ctab[i]=va*cb[i];}
    if(ctvb!=NULL) {ctvb[i]=vv*cvb;}
    /* prepare to the next loop */
    t_last=t[i]; ca_last=ca[i];
    ct1_last=ct1; ct1i_last=ct1i;
    ct2_last=ct2; ct2i_last=ct2i;
    ct3_last=ct3; ct3i_last=ct3i;
  }
 
  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;
}

simDispersion()

int simDispersion	(	double *	x,
		double *	y,
		const int	n,
		const double	tau1,
		const double	tau2,
		double *	tmp )

Simulate the effect of dispersion on a time-activity curve.

The units of rate constants must be related to the TAC time units; 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, simC1_i, corDispersion, simDelay, simTTM, simBTAC

Parameters

x	Array of sample times; must be in increasing order.
y	Array of sample values, which will be replaced here by dispersion added values.
n	Nr of samples.
tau1	First dispersion time constant (zero if no dispersion); in same time unit as sample times.
tau2	2nd dispersion time constant (zero if no dispersion); in same time unit as sample times.
tmp	Array for temporary data, for at least n samples; enter NULL to let function to allocate and free the temporary space.

Definition at line 26 of file simdispersion.c.

  {
  /* Check input */
  if(x==NULL || y==NULL || n<2) return 1;
  if(tau1<0.0 || tau2<0.0) return 2;
 
  /* Allocate memory if not allocated by user */
  double *buf;
  if(tmp!=NULL) buf=tmp; else buf=(double*)malloc(n*sizeof(double));
  if(buf==NULL) return 3;
 
  /* First dispersion */
  if(tau1>0.0) {
    double k=1.0/tau1;
    int ret=simC1(x, y, n, k, k, buf); 
    if(ret!=0) {
      if(tmp==NULL) free(buf);
      return 100+ret;
    }
    for(int i=0; i<n; i++) y[i]=buf[i];
  }
 
  /* Second dispersion */
  if(tau2>0.0) {
    double k=1.0/tau2;
    int ret=simC1(x, y, n, k, k, buf);
    if(ret!=0) {
      if(tmp==NULL) free(buf);
      return 200+ret;
    }
    for(int i=0; i<n; i++) y[i]=buf[i];
  }
 
  if(tmp==NULL) free(buf);
  return 0;
}

simIsSteadyInterval()

int simIsSteadyInterval	(	double *	x,
		const int	n,
		double *	f )

Check whether given values have steady intervals, with the first interval starting at zero. Optionally return the interval, or in case of uneven intervals, the shortest interval.

Returns

Returns 1 if values are equidistant, and 0 if intervals are uneven, or cannot be determined.

See also

convolve1D, tacInterpolateToEqualLengthFrames

Parameters

x	Array of values to check; must be sorted in ascending order.
n	Size of the data array; at least two.
f	Pointer for the interval, either the common or the shortest interval; enter NULL, if not needed.

Definition at line 66 of file convolut.c.

  {
  if(f!=NULL) *f=nan("");
  if(x==NULL || n<2) return(0);
  //if(!(fabs(x[0])>1.0E-12)) return(0); // on purpose, to catch NaN
 
  double d, freq=nan("");
  int isdif=0;
  for(int i=1; i<n; i++) {
    d=x[i]-x[i-1];
    if(!(d>1.0E-12)) return(0); // interval must be > 0
    if(isnan(freq)) {freq=d; continue;} // first interval
    if(fabs(freq-d)<1.0E-12) continue; // same interval as before
    isdif++; // different
    if(d<freq) freq=d; // save smaller interval
  }
  if(f!=NULL) *f=freq;
  if(isdif>0) return(0); // variable intervals
  /* Check that first sample time is at the midpoint of interval starting from 0 */
  d=0.5*freq; //printf("d=%g x[0]=%g\n", d, x[0]);
  if(fabs(x[0]-d)>1.0E-06) return(0);
  return(1);
}

simMBF()

int simMBF	(	double *	t,
		double *	ci,
		const int	nr,
		const double	k1,
		const double	k2,
		const double	Vfit,
		double *	ct )

Simulate myocardial tissue TAC using Iida's compartment model.

Memory for ct must be allocated in the calling program. 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, simOxygen, simSamples, simDelay

Parameters

t	Array of time values.
ci	Input activities.
nr	Number of values in TACs.
k1	Apparent k1.
k2	Apparent k2.
Vfit	Vfit
ct	Pointer for TAC array to be simulated; must be allocated.

Definition at line 28 of file sim1cm.c.

  {
  int i;
  double dt2;
  double cii, ci_last, t_last;
  double ct_last, cti, cti_last;
 
 
  /* Check for data */
  if(nr<2) return(1);
  if(t==NULL || ci==NULL || ct==NULL) return(2);
 
  /* Calculate curves */
  t_last=0.0; cii=ci_last=0.0;
  cti=ct_last=cti_last=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) {
      /* input integral */
      cii+=(ci[i]+ci_last)*dt2;
      /* Tissue compartment and its integral */
      ct[i] = (Vfit*ci[i] + k1*cii - k2*(cti_last+dt2*ct_last)) / (1.0 + dt2*k2);
      cti = cti_last + dt2*(ct_last+ct[i]);
    }
    /* set very small values to zero
    if(fabs(ct[i])<1.0e-12) ct[i]=0.0; */
    /* prepare to the next loop */
    t_last=t[i]; ci_last=ci[i];
    ct_last=ct[i]; cti_last=cti;
  }
  return(0);
}

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;
}

simRTCM()

int simRTCM	(	double *	t,
		double *	cr,
		const int	nr,
		const double	R1,
		const double	k2,
		const double	k3,
		const double	k4,
		double *	ct,
		double *	cta,
		double *	ctb )

Simulate tissue TAC using full reference tissue compartment model (original) and reference region TAC, at reference region TAC times.

Memory for ct must be allocated in the calling program. To retrieve the separate tissue compartment TACs, pointer to allocated memory for cf and/or cb can be given; if compartmental TACs are not required, NULL pointer can be given instead.

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

See also

simSRTM, simTRTM, simC2, parExampleTTACs, parExamplePerfectBolus

Returns

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

Author

Vesa Oikonen

Parameters

t	Array of time values.
cr	Reference region activities.
nr	Number of values in TACs.
R1	Ratio K1/K1'.
k2	Rate constant of the model.
k3	Rate constant of the model.
k4	Rate constant of the model.
ct	Pointer for TAC array to be simulated; must be allocated.
cta	Pointer for 1st compartment TAC to be simulated, or NULL.
ctb	Pointer for 2nd compartment TAC to be simulated, or NULL.

Definition at line 32 of file simrtcm.c.

  {
  int i;
  double f, b, w, dt2;
  double cri, cr_last, t_last;
  double cf, cf_last, cb, cb_last;
  double cfi, cfi_last, cbi, cbi_last;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(ct==NULL) return 2;
 
  /* Calculate curves */
  t_last=0.0; if(t[0]<t_last) t_last=t[0];
  cri=cr_last=0.0; cf_last=cb_last=cfi_last=cbi_last=cf=cb=cfi=cbi=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) {
      /* reference integral */
      cri+=(cr[i]+cr_last)*dt2;
      /* partial results */
      f=cfi_last+dt2*cf_last;
      b=cbi_last+dt2*cb_last;
      w=k2 + k3 + k2*k4*dt2;
      /* 1st tissue compartment and its integral */
      cf = ( (1.0 + k4*dt2)*(R1*cr[i] + k2*cri) + k4*b - w*f ) / ( 1.0 + dt2*(w+k4) );
      cfi = cfi_last + dt2*(cf_last+cf);
      /* 2nd tissue compartment and its integral */
      cb = (k3*cfi - k4*b) / (1.0 + k4*dt2);
      cbi = cbi_last + dt2*(cb_last+cb);
    }
    /* copy values to argument arrays */
    ct[i]=cf+cb;
    if(cta!=NULL) cta[i]=cf;
    if(ctb!=NULL) ctb[i]=cb;
    /* prepare to the next loop */
    t_last=t[i]; cr_last=cr[i];
    cf_last=cf; cfi_last=cfi;
    cb_last=cb; cbi_last=cbi;
  }
 
  return 0;
}

simSRTM()

int simSRTM	(	double *	t,
		double *	cr,
		const int	nr,
		const double	R1,
		const double	k2,
		const double	BP,
		double *	ct )

Simulate tissue TAC using simplified reference tissue input compartment model.

Memory for ct must be allocated in the calling program. The units of rate constants must be related to the time unit; 1/min and min, or 1/sec and sec.

See also

simRTCM, simTRTM, simC2, simC1_i

Returns

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

Author

Vesa Oikonen

Parameters

t	Array of time values.
cr	Reference region activities.
nr	Number of values in TACs.
R1	Ratio K1/K1'.
k2	Rate constant of the model.
BP	Binding potential.
ct	Pointer for TAC array to be simulated; must be allocated.

Definition at line 114 of file simrtcm.c.

  {
  int i;
  double dt2;
  double cri, cr_last, t_last;
  double ct_last, cti, cti_last;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(ct==NULL) return 2;
 
  /* Calculate curves */
  t_last=0.0; if(t[0]<t_last) t_last=t[0];
  cri=cr_last=0.0; cti=ct_last=cti_last=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) {
      /* reference integral */
      cri+=(cr[i]+cr_last)*dt2;
      /* Tissue compartment and its integral */
      ct[i] = ( R1*cr[i] + k2*cri - (k2/(1.0+BP))*(cti_last+dt2*ct_last) ) /
              ( 1.0 + dt2*(k2/(1.0+BP)) );
      cti = cti_last + dt2*(ct_last+ct[i]);
    } else { // dt==0
      ct[i]=ct_last;
    }
    /* prepare to the next loop */
    t_last=t[i]; cr_last=cr[i];
    ct_last=ct[i]; cti_last=cti;
  }
 
  return 0;
}

simTRTM()

int simTRTM	(	double *	t,
		double *	cr,
		const int	nr,
		const double	R1,
		const double	k2,
		const double	k3,
		double *	ct )

Simulate tissue TAC using reference tissue compartment model (transport limited in ref region) and reference region TAC, at reference region TAC times.

Memory for ct must be allocated in the calling program.

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

See also

simSRTM, simSRTM, simC2

Returns

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

Author

Vesa Oikonen

Parameters

t	Array of time values.
cr	Reference region activities.
nr	Number of values in TACs.
R1	Ratio K1/K1'.
k2	Rate constant of the model.
k3	Rate constant of the model.
ct	Pointer for TAC array to be simulated; must be allocated.

Definition at line 182 of file simrtcm.c.

  {
  int i;
  double dt2;
  double cri, cr_last, t_last;
  double ct_last, cti, cti_last;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(ct==NULL) return 2;
 
  /* Calculate curves */
  t_last=0.0; if(t[0]<t_last) t_last=t[0];
  cri=cr_last=0.0; cti=ct_last=cti_last=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) {
      /* reference integral */
      cri+=(cr[i]+cr_last)*dt2;
      /* Tissue compartment and its integral */
      ct[i] = ( R1*cr[i] + R1*k3*cri - (k2+k3)*(cti_last+dt2*ct_last) ) / ( 1.0 + dt2*(k2+k3) );
      cti = cti_last + dt2*(ct_last+ct[i]);
    } else { // dt==0
      ct[i]=ct_last;
    }
    /* prepare to the next loop */
    t_last=t[i]; cr_last=cr[i];
    ct_last=ct[i]; cti_last=cti;
  }
 
  return 0;
}

simTTM()

int simTTM	(	double *	t,
		double *	c0,
		const int	n,
		const double	k,
		const int	cn,
		double *	cout )

Simulate output of n-compartmental transit-time model, at input TAC sample times.

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

See also

simTTM_i, simDispersion, simC1, simSamples, simDelay

Returns

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

Author

Vesa Oikonen

Parameters

t	Array of time values.
c0	Array of input values.
n	Number of values in TACs.
k	Rate constant of the model.
cn	Number of compartments (cn>0).
cout	Pointer for TAC array to be simulated; must be allocated.

Definition at line 144 of file simdispersion.c.

  {
  /* Check for data */
  if(n<2 || cn<1) return(1);
  if(t==NULL || c0==NULL || cout==NULL) return(2);
  /* Check the rate constant */
  if(!(k>=0.0)) return(3);
 
  double c[cn+1], ci[cn+1]; // Compartmental concentrations for current sample
  double c_last[cn+1], ci_last[cn+1]; // Compartmental concentrations for previous sample
  for(int j=0; j<=cn; j++) c[j]=c_last[j]=ci[j]=ci_last[j]=0.0;
 
  /* Calculate curves */
#if(0)
  printf("\nTime\tC0\tiC0");
  for(int j=1; j<=cn; j++) printf("\tC%d", j);
  printf("\n");
#endif
  double t_last=0.0; if(t[0]<t_last) t_last=t[0];
  double c0_last=0.0, c0i=0.0;
  for(int i=0; i<n; i++) { // loop through sample times
    /* delta time / 2 */
    double dt2=0.5*(t[i]-t_last); if(!(dt2>=0.0)) return(5);
    /* input integral */
    c0i+=(c0[i]+c0_last)*dt2;
    /* k/(1+k*(dt/2)) */
    double kdt=k/(1.0+dt2*k);
    /* calculate compartmental concentrations */
    ci[0]=c0i; 
    if(dt2>0.0) {
      for(int j=1; j<=cn; j++) {
        c[j] = kdt * (ci[j-1] - ci_last[j] - dt2*c_last[j]);
        ci[j] = ci_last[j] + dt2*(c_last[j]+c[j]);
      }
    } else { // sample time same as previously, thus concentration do not change either
      for(int j=0; j<=cn; j++) {
        c[j]=c_last[j];
        ci[j]=ci_last[j];
      }
    }
#if(0)
    printf("%g\t%g\t%g", t[i], c0[i], ci[0]);
    for(int j=1; j<=cn; j++) printf("\t%g", c[j]);
    printf("\n");
#endif
    cout[i]=c[cn];
    /* prepare to the next loop */
    t_last=t[i]; c0_last=c0[i];
    for(int j=0; j<=cn; j++) {
      c_last[j]=c[j];
      ci_last[j]=ci[j];
    }
  }
 
  return(0);
}

simTTM_i()

int simTTM_i	(	double *	t,
		double *	c0i,
		const int	n,
		const double	k,
		const int	cn,
		double *	cout )

Simulate output of n-compartmental transit-time model, at input TAC sample times.

This version uses integral of input function, enabling user to fully control the calculation of integral, which is more precise when input is based on an integrable mathematical function.

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

See also

simTTM, simDispersion, simC1_i, simSamples

Returns

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

Author

Vesa Oikonen

Parameters

t	Array of time values.
c0i	Array of AUC 0-t of input function.
n	Number of values in TACs.
k	Rate constant of the model.
cn	Number of compartments (cn>0).
cout	Pointer for TAC array to be simulated; must be allocated.

Definition at line 227 of file simdispersion.c.

  {
  /* Check for data */
  if(n<2 || cn<1) return(1);
  if(t==NULL || c0i==NULL || cout==NULL) return(2);
  /* Check the rate constant */
  if(!(k>=0.0)) return(3);
 
  double c[cn+1], ci[cn+1]; // Compartmental concentrations for current sample
  double c_last[cn+1], ci_last[cn+1]; // Compartmental concentrations for previous sample
  for(int j=0; j<=cn; j++) c[j]=c_last[j]=ci[j]=ci_last[j]=0.0;
 
  /* Calculate curves */
#if(0)
  printf("\nTime\tiC0");
  for(int j=1; j<=cn; j++) printf("\tC%d", j);
  printf("\n");
#endif
  double t_last=0.0; if(t[0]<t_last) t_last=t[0];
  for(int i=0; i<n; i++) { // loop through sample times
    /* delta time / 2 */
    double dt2=0.5*(t[i]-t_last); if(!(dt2>=0.0)) return(5);
    /* k/(1+k*(dt/2)) */
    double kdt=k/(1.0+dt2*k);
    /* calculate compartmental concentrations */
    ci[0]=c0i[i]; 
    if(dt2>0.0) {
      for(int j=1; j<=cn; j++) {
        c[j] = kdt * (ci[j-1] - ci_last[j] - dt2*c_last[j]);
        ci[j] = ci_last[j] + dt2*(c_last[j]+c[j]);
      }
    } else { // sample time same as previously, thus concentration do not change either
      for(int j=0; j<=cn; j++) {
        c[j]=c_last[j];
        ci[j]=ci_last[j];
      }
    }
#if(0)
    printf("%g\t%g", t[i], ci[0]);
    for(int j=1; j<=cn; j++) printf("\t%g", c[j]);
    printf("\n");
#endif
    cout[i]=c[cn];
    /* prepare to the next loop */
    t_last=t[i];
    for(int j=0; j<=cn; j++) {
      c_last[j]=c[j];
      ci_last[j]=ci[j];
    }
  }
 
  return(0);
}