simblood.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
/*****************************************************************************/
#include "tpccm.h"
/*****************************************************************************/
 
/*****************************************************************************/
void icmparcInit(
  ICMPARC *d
) {
  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;
}
/*****************************************************************************/
 
/*****************************************************************************/
int icmparcAddMetabolites(
  ICMPARC *d,
  unsigned const int mNr
) {
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/
int icmparcAllocateTACs(
  ICMPARC *d,
  unsigned const int sNr,
  const int sub
) {
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/
void icmparcFree(
  ICMPARC *d
) {
  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;
}
/*****************************************************************************/
 
/*****************************************************************************/
int simBTAC(
  double *t,
  const unsigned int nr,
  ICMPARC *p,
  double *c_BA
) {
  /* 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);
}
/*****************************************************************************/
 
/*****************************************************************************/
int simBM1(
  double *t,
  double *ct,
  const int nr,
  const double k,
  double *cp
) {
  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;
}
/*****************************************************************************/
 
/*****************************************************************************/