TPCCLIB
Loading...
Searching...
No Matches
simdicm.c
Go to the documentation of this file.
1
4/*****************************************************************************/
5#include "tpcclibConfig.h"
6/*****************************************************************************/
7#include <stdio.h>
8#include <stdlib.h>
9#include <math.h>
10#include <time.h>
11#include <string.h>
12/*****************************************************************************/
13#include "tpccm.h"
14/*****************************************************************************/
15
16/*****************************************************************************/
30int simC1DI(
32 double *t,
34 double *cba,
36 double *cbb,
38 const int nr,
40 const double k1a,
42 const double k1b,
44 const double k2,
46 double *ct
47) {
48 /* Check for data */
49 if(nr<2 || t==NULL || cba==NULL || cbb==NULL) return(1);
50 if(ct==NULL) return(2);
51
52 /* Check parameters */
53 if(!(k1a>=0.0)) return(3);
54 if(!(k1b>=0.0)) return(4);
55 if(!(k2>=0.0)) return(5);
56
57 /* Calculate curves */
58 double t_last=0.0; if(t[0]<t_last) t_last=t[0];
59 double cbai=0.0, cba_last=0.0, cbbi=0.0, cbb_last=0.0;
60 double ct1=0.0, ct1i=0.0, ct1_last=0.0, ct1i_last=0.0;
61 for(int i=0; i<nr; i++) {
62 /* delta time / 2 */
63 double dt2=0.5*(t[i]-t_last);
64 /* calculate values */
65 if(dt2<0.0) {
66 return(6);
67 } else if(dt2>0.0) {
68 /* input integrals */
69 cbai+=(cba[i]+cba_last)*dt2;
70 cbbi+=(cbb[i]+cbb_last)*dt2;
71 /* tissue compartment and its integral */
72 ct1 = (k1a*cbai + k1b*cbbi - k2*(ct1i_last+dt2*ct1_last)) / (1.0 + dt2*k2);
73 ct1i = ct1i_last + dt2*(ct1_last+ct1);
74 }
75 /* copy values to argument array */
76 ct[i]=ct1;
77 /* prepare to the next loop */
78 t_last=t[i]; cba_last=cba[i]; cbb_last=cbb[i];
79 ct1_last=ct1; ct1i_last=ct1i;
80 }
81
82 return(0);
83}
84/*****************************************************************************/
85
86/*****************************************************************************/
101int simC3DIvs(
103 double *t,
105 double *ca1,
107 double *ca2,
109 double *cb,
111 const int nr,
113 const double k1,
115 const double k2,
117 const double k3,
119 const double k4,
121 const double k5,
123 const double k6,
125 const double k1b,
127 const double k2b,
129 const double f,
131 const double vb,
133 const double fa,
137 const int vvm,
139 double *scpet,
141 double *sct1,
143 double *sct2,
145 double *sct3,
147 double *sct1b,
149 double *sctab,
151 double *sctvb
152) {
153 int i;
154 double b, c, d, e, w, z, dt2, va, vv;
155 double ca1i, ca1_last, ca2i, ca2_last, t_last, dct, cvb;
156 double ct1, ct1_last, ct2, ct2_last, ct3, ct3_last;
157 double ct1i, ct1i_last, ct2i, ct2i_last, ct3i, ct3i_last;
158 double ct1b, ct1b_last, ct1bi, ct1bi_last;
159
160
161 /* Check for data */
162 if(nr<2) return 1;
163 if(scpet==NULL) return 2;
164
165 /* Check parameters */
166 if(k1<0.0) return 3;
167 if(vb<0.0 || vb>=1.0) return 4;
168 if(fa<=0.0 || fa>1.0) return 5;
169 va=fa*vb; vv=(1.0-fa)*vb;
170
171 /* Calculate curves */
172 t_last=0.0; if(t[0]<t_last) t_last=t[0];
173 ca1i=ca1_last=ca2i=ca2_last=0.0;
174 ct1_last=ct2_last=ct3_last=ct1i_last=ct2i_last=ct3i_last=0.0;
175 ct1b_last=ct1bi_last=0.0;
176 ct1=ct2=ct3=ct1b=ct1i=ct2i=ct3i=ct1bi=0.0;
177 for(i=0; i<nr; i++) {
178 /* delta time / 2 */
179 dt2=0.5*(t[i]-t_last);
180 /* calculate values */
181 if(dt2<0.0) {
182 return 5;
183 } else if(dt2>0.0) {
184 /* arterial integrals */
185 ca1i+=(ca1[i]+ca1_last)*dt2;
186 ca2i+=(ca2[i]+ca2_last)*dt2;
187 /* partial results */
188 b=ct1i_last+dt2*ct1_last;
189 c=ct2i_last+dt2*ct2_last;
190 d=ct3i_last+dt2*ct3_last;
191 e=ct1bi_last+dt2*ct1b_last;
192 w=k4 + k5 - (k5*k6*dt2)/(1.0+k6*dt2);
193 z=1.0+w*dt2;
194 /* 1st tissue compartment and its integral */
195 ct1 = (
196 + k1*z*ca1i + (k3*k4*dt2 - (k2+k3)*z)*b
197 + k4*c + k4*k6*dt2*d/(1.0+k6*dt2)
198 ) / ( z*(1.0 + dt2*(k2+k3)) - k3*k4*dt2*dt2 );
199 ct1i = ct1i_last + dt2*(ct1_last+ct1);
200 ct1b = (k1b*ca2i - k2b*e) / (1.0 + dt2*k2b);
201 ct1bi = ct1bi_last + dt2*(ct1b_last+ct1b);
202 /* 2nd tissue compartment and its integral */
203 ct2 = (k3*ct1i - w*c + k6*d/(1.0+k6*dt2)) / z;
204 ct2i = ct2i_last + dt2*(ct2_last+ct2);
205 /* 3rd tissue compartment and its integral */
206 ct3 = (k5*ct2i - k6*d) / (1.0 + k6*dt2);
207 ct3i = ct3i_last + dt2*(ct3_last+ct3);
208 }
209 /* Venous curve */
210 if(f>0.) {
211 dct = k1*ca1[i] - k2*ct1 + k1b*ca2[i] - k2b*ct1b;
212 cvb = cb[i] - dct/f;
213 } else cvb=cb[i];
214 /* copy values to argument arrays */
215 if(vvm==0) scpet[i]= va*cb[i] + vv*cvb + (1.0-vb)*(ct1+ct2+ct3+ct1b);
216 else scpet[i]= va*cb[i] + vv*cvb + (ct1+ct2+ct3+ct1b);
217 if(sct1!=NULL) {sct1[i]=(1.0-vb)*ct1; if(vvm==0) sct1[i]*=(1.0-vb);}
218 if(sct2!=NULL) {sct2[i]=(1.0-vb)*ct2; if(vvm==0) sct2[i]*=(1.0-vb);}
219 if(sct3!=NULL) {sct3[i]=(1.0-vb)*ct3; if(vvm==0) sct3[i]*=(1.0-vb);}
220 if(sct1b!=NULL) {sct1b[i]=(1.0-vb)*ct1b; if(vvm==0) sct1b[i]*=(1.0-vb);}
221 if(sctab!=NULL) sctab[i]=va*cb[i];
222 if(sctvb!=NULL) sctvb[i]=vv*cvb;
223 /* prepare to the next loop */
224 t_last=t[i]; ca1_last=ca1[i]; ca2_last=ca2[i];
225 ct1_last=ct1; ct1i_last=ct1i;
226 ct2_last=ct2; ct2i_last=ct2i;
227 ct3_last=ct3; ct3i_last=ct3i;
228 ct1b_last=ct1b; ct1bi_last=ct1bi;
229 }
230
231 return 0;
232}
233/*****************************************************************************/
234
235/*****************************************************************************/
251int simC4DIvp(
253 double *t,
255 double *ca1,
257 double *ca2,
259 double *cb,
261 const int nr,
263 const double k1,
265 const double k2,
267 const double k3,
269 const double k4,
271 const double k5,
273 const double k6,
275 const double k7,
277 const double km,
279 double k1b,
281 const double k2b,
283 const double f,
285 const double vb,
287 double fa,
291 const int vvm,
293 double *scpet,
295 double *sct1,
297 double *sct2,
299 double *sct3,
301 double *sct1b,
303 double *sctab,
305 double *sctvb,
307 const int verbose
308) {
309 int i;
310 double b, c, d, e, pt, qt, dt2, va, vv;
311 double ca1i, ca1_last, ca2i, ca2_last, t_last, dct, cvb;
312 double ct1, ct1_last, ct2, ct2_last, ct3, ct3_last;
313 double ct1i, ct1i_last, ct2i, ct2i_last, ct3i, ct3i_last;
314 double ct1b, ct1b_last, ct1bi, ct1bi_last;
315
316
317 if(verbose>0) {
318 printf("%s()\n", __func__);
319 if(verbose>1) {
320 printf(" k1 := %g\n", k1);
321 printf(" k2 := %g\n", k2);
322 printf(" k3 := %g\n", k3);
323 printf(" k4 := %g\n", k4);
324 printf(" k5 := %g\n", k5);
325 printf(" k6 := %g\n", k6);
326 printf(" k7 := %g\n", k7);
327 printf(" km := %g\n", km);
328 printf(" k1b := %g\n", k1b);
329 printf(" k2b := %g\n", k2b);
330 printf(" vb := %g\n", vb);
331 printf(" fa := %g\n", fa);
332 printf(" f := %g\n", f);
333 }
334 }
335
336 /* Check for data */
337 if(nr<2) return 1;
338 if(scpet==NULL) return 2;
339
340 /* Check parameters */
341 if(k1<0.0 || k1b<0.0) return 3;
342 if(vb<0.0 || vb>=1.0) return 4;
343 if(fa<0.0 || fa>1.0) return 5;
344 va=fa*vb; vv=(1.0-fa)*vb;
345
346 /* Calculate curves */
347 t_last=0.0; if(t[0]<t_last) t_last=t[0];
348 ca1i=ca1_last=ca2i=ca2_last=0.0;
349 ct1_last=ct2_last=ct3_last=ct1i_last=ct2i_last=ct3i_last=ct1b_last=
350 ct1bi_last=0.0;
351 ct1=ct2=ct3=ct1b=ct1i=ct2i=ct3i=ct1bi=0.0;
352 for(i=0; i<nr; i++) {
353 /* delta time / 2 */
354 dt2=0.5*(t[i]-t_last);
355 /* calculate values */
356 if(dt2<0.0) {
357 return 5;
358 } else if(dt2>0.0) {
359 /* arterial integrals */
360 ca1i+=(ca1[i]+ca1_last)*dt2;
361 ca2i+=(ca2[i]+ca2_last)*dt2;
362 /* partial results */
363 b=ct1i_last+dt2*ct1_last;
364 c=ct2i_last+dt2*ct2_last;
365 d=ct3i_last+dt2*ct3_last;
366 e=ct1bi_last+dt2*ct1b_last;
367 pt=k6+k7;
368 qt=k2+k3+k5+km-(k3*k4*dt2)/(1.0+k4*dt2)-(k5*k6*dt2)/(1.0+pt*dt2);
369 /* 1st tissue compartment and its integral */
370 ct1 = (k1/(1.0+qt*dt2))*ca1i
371 - (qt/(1.0+qt*dt2))*b
372 + (k4/((1.0+qt*dt2)*(1.0+k4*dt2)))*c
373 + (k6/((1.0+qt*dt2)*(1.0+pt*dt2)))*d;
374 ct1i = ct1i_last + dt2*(ct1_last+ct1);
375 /* 2nd tissue compartment and its integral */
376 ct2 = (k3/(1.0+k4*dt2))*ct1i
377 - (k4/(1.0+k4*dt2))*c;
378 ct2i = ct2i_last + dt2*(ct2_last+ct2);
379 /* 3rd tissue compartment and its integral */
380 ct3 = (k5/(1.0+pt*dt2))*ct1i
381 - (pt/(1.0+pt*dt2))*d;
382 ct3i = ct3i_last + dt2*(ct3_last+ct3);
383 /* 4th tissue compartment (the 1st for tracer 2) and its integral */
384 ct1b = (k1b/(1.0+k2b*dt2))*ca2i
385 - (k2b/(1.0+k2b*dt2))*e
386 + (km/(1.0+k2b*dt2))*ct1i;
387 ct1bi = ct1bi_last + dt2*(ct1b_last+ct1b);
388 }
389 /* Venous curve */
390 if(f>0.) {
391 dct = k1*ca1[i] - k2*ct1 - k7*ct3 + k1b*ca2[i] - k2b*ct1b;
392 cvb = cb[i] - dct/f;
393 } else cvb=cb[i];
394 /* copy values to argument arrays */
395 if(vvm==0) scpet[i]= va*cb[i] + vv*cvb + (1.0-vb)*(ct1+ct2+ct3+ct1b);
396 else scpet[i]= va*cb[i] + vv*cvb + (ct1+ct2+ct3+ct1b);
397 if(sct1!=NULL) {sct1[i]=(1.0-vb)*ct1; if(vvm==0) sct1[i]*=(1.0-vb);}
398 if(sct2!=NULL) {sct2[i]=(1.0-vb)*ct2; if(vvm==0) sct2[i]*=(1.0-vb);}
399 if(sct3!=NULL) {sct3[i]=(1.0-vb)*ct3; if(vvm==0) sct3[i]*=(1.0-vb);}
400 if(sct1b!=NULL) {sct1b[i]=(1.0-vb)*ct1b; if(vvm==0) sct1b[i]*=(1.0-vb);}
401 if(sctab!=NULL) sctab[i]=va*cb[i];
402 if(sctvb!=NULL) sctvb[i]=vv*cvb;
403 /* prepare to the next loop */
404 t_last=t[i]; ca1_last=ca1[i]; ca2_last=ca2[i];
405 ct1_last=ct1; ct1i_last=ct1i;
406 ct2_last=ct2; ct2i_last=ct2i;
407 ct3_last=ct3; ct3i_last=ct3i;
408 ct1b_last=ct1b; ct1bi_last=ct1bi;
409 }
410
411 if(verbose>2) {
412 printf("AUC 0-%g:\n", t_last);
413 printf(" ca1i := %g\n", ca1i);
414 printf(" ca2i := %g\n", ca2i);
415 printf(" ct1i := %g\n", ct1i_last);
416 printf(" ct2i := %g\n", ct2i_last);
417 printf(" ct3i := %g\n", ct3i_last);
418 printf(" ct1bi := %g\n", ct1bi_last);
419 }
420
421 return 0;
422}
423/*****************************************************************************/
424
425/*****************************************************************************/
441int simC4DIvs(
443 double *t,
445 double *ca1,
447 double *ca2,
449 double *cb,
451 const int nr,
453 const double k1,
455 const double k2,
457 const double k3,
459 const double k4,
461 const double k5,
463 const double k6,
465 const double k7,
467 const double km,
469 double k1b,
471 double k2b,
473 double f,
475 double vb,
477 double fa,
481 const int vvm,
483 double *scpet,
485 double *sct1,
487 double *sct2,
489 double *sct3,
491 double *sct1b,
493 double *sctab,
495 double *sctvb,
497 const int verbose
498) {
499 int i;
500 double b, c, d, e, pt, qt, rt, dt2, va, vv;
501 double ca1i, ca1_last, ca2i, ca2_last, t_last, dct, cvb;
502 double ct1, ct1_last, ct2, ct2_last, ct3, ct3_last;
503 double ct1i, ct1i_last, ct2i, ct2i_last, ct3i, ct3i_last;
504 double ct1b, ct1b_last, ct1bi, ct1bi_last;
505
506
507 if(verbose>0) {
508 printf("%s()\n", __func__);
509 if(verbose>1) {
510 printf(" k1 := %g\n", k1);
511 printf(" k2 := %g\n", k2);
512 printf(" k3 := %g\n", k3);
513 printf(" k4 := %g\n", k4);
514 printf(" k5 := %g\n", k5);
515 printf(" k6 := %g\n", k6);
516 printf(" k7 := %g\n", k7);
517 printf(" km := %g\n", km);
518 printf(" k1b := %g\n", k1b);
519 printf(" k2b := %g\n", k2b);
520 printf(" vb := %g\n", vb);
521 printf(" fa := %g\n", fa);
522 printf(" f := %g\n", f);
523 }
524 }
525
526 /* Check for data */
527 if(nr<2) return 1;
528 if(scpet==NULL) return 2;
529
530 /* Check parameters */
531 if(k1<0.0 || k1b<0.0) return 3;
532 if(vb<0.0 || vb>=1.0) return 4;
533 if(fa<0.0 || fa>1.0) return 5;
534 va=fa*vb; vv=(1.0-fa)*vb;
535
536 /* Calculate curves */
537 t_last=0.0; if(t[0]<t_last) t_last=t[0];
538 ca1i=ca1_last=ca2i=ca2_last=0.0;
539 ct1_last=ct2_last=ct3_last=ct1i_last=ct2i_last=ct3i_last=ct1b_last=
540 ct1bi_last=0.0;
541 ct1=ct2=ct3=ct1b=ct1i=ct2i=ct3i=ct1bi=0.0;
542 for(i=0; i<nr; i++) {
543 /* delta time / 2 */
544 dt2=0.5*(t[i]-t_last);
545 /* calculate values */
546 if(dt2<0.0) {
547 return 5;
548 } else if(dt2>0.0) {
549 /* arterial integrals */
550 ca1i+=(ca1[i]+ca1_last)*dt2;
551 ca2i+=(ca2[i]+ca2_last)*dt2;
552 //printf("ca1[%d]=%g int=%g\n", i, ca1[i], ca1i);
553 //printf("ca2[%d]=%g int=%g\n", i, ca2[i], ca2i);
554 /* partial results */
555 b=ct1i_last+dt2*ct1_last;
556 c=ct2i_last+dt2*ct2_last;
557 d=ct3i_last+dt2*ct3_last;
558 e=ct1bi_last+dt2*ct1b_last;
559 pt=k6+k7;
560 qt=k4+k5-(k5*k6*dt2)/(1.0+pt*dt2);
561 rt=k2+k3+km-(k3*k4*dt2)/(1.0+qt*dt2);
562 /* 1st tissue compartment and its integral */
563 ct1 = (k1/(1.0+rt*dt2))*ca1i
564 - (rt/(1.0+rt*dt2))*b
565 + (k4/((1.0+qt*dt2)*(1.0+rt*dt2)))*c
566 + ((k4*k6*dt2)/((1.0+pt*dt2)*(1.0+qt*dt2)*(1.0+rt*dt2)))*d;
567 ct1i = ct1i_last + dt2*(ct1_last+ct1);
568 /* 2nd tissue compartment and its integral */
569 ct2 = (k3/(1.0+qt*dt2))*ct1i
570 - (qt/(1.0+qt*dt2))*c
571 + (k6/((1.0+pt*dt2)*(1.0+qt*dt2)))*d;
572 ct2i = ct2i_last + dt2*(ct2_last+ct2);
573 /* 3rd tissue compartment and its integral */
574 ct3 = (k5/(1.0+pt*dt2))*ct2i
575 - (pt/(1.0+pt*dt2))*d;
576 ct3i = ct3i_last + dt2*(ct3_last+ct3);
577 /* 4th tissue compartment (the 1st for tracer 2) and its integral */
578 ct1b = (k1b/(1.0+k2b*dt2))*ca2i
579 - (k2b/(1.0+k2b*dt2))*e
580 + (km/(1.0+k2b*dt2))*ct1i;
581 ct1bi = ct1bi_last + dt2*(ct1b_last+ct1b);
582 }
583 /* Venous curve */
584 if(f>0.) {
585 dct = k1*ca1[i] - k2*ct1 - k7*ct3 + k1b*ca2[i] - k2b*ct1b;
586 cvb = cb[i] - dct/f;
587 } else cvb=cb[i];
588 /* copy values to argument arrays; set very small values to zero */
589 if(vvm==0) scpet[i]= va*cb[i] + vv*cvb + (1.0-vb)*(ct1+ct2+ct3+ct1b);
590 else scpet[i]= va*cb[i] + vv*cvb + (ct1+ct2+ct3+ct1b);
591 if(sct1!=NULL) {sct1[i]=(1.0-vb)*ct1; if(vvm==0) sct1[i]*=(1.0-vb);}
592 if(sct2!=NULL) {sct2[i]=(1.0-vb)*ct2; if(vvm==0) sct2[i]*=(1.0-vb);}
593 if(sct3!=NULL) {sct3[i]=(1.0-vb)*ct3; if(vvm==0) sct3[i]*=(1.0-vb);}
594 if(sct1b!=NULL) {sct1b[i]=(1.0-vb)*ct1b; if(vvm==0) sct1b[i]*=(1.0-vb);}
595 if(sctab!=NULL) sctab[i]=va*cb[i];
596 if(sctvb!=NULL) sctvb[i]=vv*cvb;
597 /* prepare to the next loop */
598 t_last=t[i]; ca1_last=ca1[i]; ca2_last=ca2[i];
599 ct1_last=ct1; ct1i_last=ct1i;
600 ct2_last=ct2; ct2i_last=ct2i;
601 ct3_last=ct3; ct3i_last=ct3i;
602 ct1b_last=ct1b; ct1bi_last=ct1bi;
603 }
604
605 if(verbose>2) {
606 printf("AUC 0-%g:\n", t_last);
607 printf(" ca1i := %g\n", ca1i);
608 printf(" ca2i := %g\n", ca2i);
609 printf(" ct1i := %g\n", ct1i_last);
610 printf(" ct2i := %g\n", ct2i_last);
611 printf(" ct3i := %g\n", ct3i_last);
612 printf(" ct1bi := %g\n", ct1bi_last);
613 }
614
615 return 0;
616}
617/*****************************************************************************/
618
619/*****************************************************************************/
simC1DI
int simC1DI(double *t, double *cba, double *cbb, const int nr, const double k1a, const double k1b, const double k2, double *ct)
Definition simdicm.c:30
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)
Definition simdicm.c:251
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)
Definition simdicm.c:101
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)
Definition simdicm.c:441
tpccm.h
Header file for libtpccm.
Generated on for TPCCLIB by doxygen 1.15.0