Help on a cell transit model

From: "Bonate, Peter" Peter.Bonate@genzyme.com Subject: [NMusers] Help on a cell transit model Date: Tue, 28 Feb 2006 11:08:29 -0500 Dear all, I was wondering if someone could help me with a problem I am having in setting up Friberg's cell transit model of cell count kinetics. Here is the control stream I am using. For simplicity I removed the feedback part to the model. The PK data is fixed to the best 2-compartment model. $PROBLEM CELL TRANSIT MODEL $DATA TRANSIT_RED.DAT IGNORE=C WIDE $INPUT ID DATE=DROP TIME DV2 AMT RATE MDV EVID CMT TRT VIST DAY X2 PL X1 SEX DV $SUBROUTINE ADVAN9 TOL=5 $MODEL NCOMP=7 COMP=(CENTRAL) COMP=(PERIP) COMP=(STEM) COMP=(TRANSIT1) COMP=(TRANSIT2) COMP=(TRANSIT3) COMP=(CELL) $PK CL = THETA(1)*(X1/30)**0.75*(X2/0.6)**THETA(5) TVV1 = THETA(2) TVQ2 = THETA(3) TVV2 = THETA(4) IF (VIST .EQ. 0) THEN V1 = TVV1*EXP(ETA(1) + ETA(2)) Q = TVQ2*EXP(ETA(5)) V2 = TVV2 END IF IF (VIST .EQ. 12) THEN V1 = TVV1*THETA(6)*EXP(ETA(1) + ETA(3)) Q = TVQ2*EXP(ETA(6)) V2 = TVV2 END IF IF (VIST .EQ. 24) THEN V1 = TVV1*THETA(6)*EXP(ETA(1) + ETA(4)) Q = TVQ2*EXP(ETA(7)) V2 = TVV2 END IF S1 = V1 BASE = THETA(7)*EXP(ETA(8)) MTT = THETA(8)*EXP(ETA(9)) KP = 4/MTT F3 = BASE F4 = BASE F5 = BASE F6 = BASE F7 = BASE SLOP = THETA(9)*EXP(ETA(10)) $ERROR PKY=0 PDY=0 IPRED = F IF (F .GT. 0 .AND. CMT .EQ. 1 .AND. TIME .GT. 0) IPRED = LOG(F) IF (CMT .EQ. 1 .AND. TIME .GT. 0) PKY = IPRED + EPS(1) IF (CMT .EQ. 3 .AND. TIME .GT. 0) PDY = IPRED*EXP(EPS(2)) + EPS(3) Y = PKY + PDY $DES CP = A(1)/V1 DRUG = SLOP*CP DADT(1) = -(CL + Q)*A(1)/V1 + Q*A(2)/V2 DADT(2) = Q*A(1)/V1 - Q*A(2)/V2 DADT(3) = KP*A(3) - KP*A(3) ; STEM DADT(4) = KP*A(3) - KP*A(4) ; TRANSIT 1 DADT(5) = KP*A(4) - KP*A(5) ; TRANSIT 2 DADT(6) = KP*A(5) - KP*A(6) ; TRANSIT 3 DADT(7) = KP*A(6) - KP*(1 + DRUG)*A(7) ; CELL COUNTS A3 = A(3) A4 = A(4) A5 = A(5) A6 = A(6) A7 = A(7) $THETA (16.5 FIXED); CLEARANCE (ML/H) (6560 FIXED); CENTRAL VOLUME (ML) (14.1 FIXED); INTERCOMPARTMENTAL CLEARANCE (8640 FIXED); PERIPHERAL VOLUME (ML) (-0.357 FIXED); EFFECT OF X2 ON CL (0.558 FIXED); EFFECT OF VISIT ON V1 (0, 1.6); BASELINE CELL COUNT (0, 200); MEAN TRANSIT TIME (0, 1); SLOPE EFFECT FOR CONCENTRATION $OMEGA 0.0491 FIXED $OMEGA BLOCK(1) 0.0585 FIXED $OMEGA BLOCK(1) SAME $OMEGA BLOCK(1) SAME $OMEGA BLOCK(1) 0.0723 FIXED $OMEGA BLOCK(1) SAME $OMEGA BLOCK(1) SAME $OMEGA 0.1 0.1 0.1 $SIGMA 0.0539 FIXED $SIGMA 0.1 0.01 ;$EST NOABORT METHOD=CONDITIONAL SIGDIGITS=3 PRINT=10 POSTHOC MAXEVAL=9999 ;$COV $SIMULATION (23523)(4466) ONLYSIMULATION $TABLE FILE=PD2.TXT NOPRINT NOHEADER ID TIME DV AMT RATE MDV EVID CMT TRT VIST DAY X2 PL X1 CL V1 Q V2 BASE MTT SLOP IPRED CP ETA(8) ETA(9) ETA(10) DRUG A3 A4 A5 A6 A7 Here are the first 30 records of the 1st individual: Obs id siteptid visit Day date time dt evid mdv dv cmt amt rate X2 X1 1 1004 201-1004 0 1 01/20/2003 11:15 1.000 1 1 . 1 12000 3000 1.72 71.34 2 1004 201-1004 0 1 01/20/2003 11:15 1.000 1 1 . 3 1 0 1.72 71.34 3 1004 201-1004 0 1 01/20/2003 11:15 1.000 1 1 . 4 1 0 1.72 71.34 4 1004 201-1004 0 1 01/20/2003 11:15 1.000 1 1 . 5 1 0 1.72 71.34 5 1004 201-1004 0 1 01/20/2003 11:15 1.000 1 1 . 6 1 0 1.72 71.34 6 1004 201-1004 0 1 01/20/2003 11:15 1.000 1 1 . 7 1 0 1.72 71.34 7 1004 201-1004 0 1 01/20/2003 15:15 1.167 0 0 -0.47804 1 . . 1.72 71.34 8 1004 201-1004 0 1 01/20/2003 18:00 1.281 0 0 1.72000 7 . . 1.72 71.34 9 1004 201-1004 0 1 01/20/2003 19:15 1.333 0 0 -0.54473 1 . . 1.72 71.34 10 1004 201-1004 0 2 01/21/2003 10:14 1.958 1 1 . 1 12000 3000 1.51 71.39 11 1004 201-1004 0 2 01/21/2003 14:14 2.124 0 0 0.37844 1 . . 1.51 71.39 12 1004 201-1004 0 2 01/21/2003 18:14 2.291 0 0 0.56531 1 . . 1.51 71.39 13 1004 201-1004 0 3 01/22/2003 10:11 2.956 0 0 0.26236 1 . . 1.30 71.44 14 1004 201-1004 0 3 01/22/2003 10:14 2.958 1 1 . 1 12000 3000 1.30 71.44 15 1004 201-1004 0 3 01/22/2003 14:14 3.124 0 0 1.49065 1 . . 1.30 71.44 16 1004 201-1004 0 3 01/22/2003 18:18 3.294 0 0 1.54330 1 . . 1.30 71.44 17 1004 201-1004 0 4 01/23/2003 9:46 3.938 0 0 1.30833 1 . . 1.09 71.49 18 1004 201-1004 0 4 01/23/2003 10:23 3.964 1 1 . 1 12000 3934 1.09 71.49 19 1004 201-1004 0 4 01/23/2003 14:24 4.131 0 0 1.67147 1 . . 1.09 71.49 20 1004 201-1004 0 4 01/23/2003 18:25 4.299 0 0 1.55181 1 . . 1.09 71.49 21 1004 201-1004 0 5 01/24/2003 9:40 4.934 0 0 1.36098 1 . . 0.88 71.54 22 1004 201-1004 0 5 01/24/2003 10:17 4.960 1 1 . 1 12000 6000 0.88 71.54 23 1004 201-1004 0 5 01/24/2003 14:15 5.125 0 0 1.65632 1 . . 0.88 71.54 24 1004 201-1004 0 5 01/24/2003 18:14 5.291 0 0 1.69745 1 . . 0.88 71.54 25 1004 201-1004 0 9 01/28/2003 15:15 9.167 0 0 1.24703 1 . . 0.03 71.74 26 1004 201-1004 0 9 01/28/2003 18:00 9.281 0 0 0.03000 7 . . 0.03 71.74 27 1004 201-1004 0 16 02/04/2003 15:15 16.167 0 0 -0.06188 1 . . 0.05 72.09 28 1004 201-1004 0 30 02/18/2003 18:00 30.281 0 0 0.10000 7 . . 0.10 72.80 29 1004 201-1004 0 79 04/08/2003 18:00 79.281 0 0 0.37000 7 . . 0.37 72.80 30 1004 201-1004 0 177 07/15/2003 18:00 177.281 0 0 0.46000 7 . . 0.46 72.80 When I run the control stream in simulation mode, the very first record at time 0 for the first subject gives initial estimates in all cell count compartments (A3 to A7) as zero. Also, the next subject (1006) has at time zero, not 0 amounts in the compartments, but the last values for the subject prior. Obs ID visit days TIME IPRED DV EVID CMT MDV A3 A4 A5 A6 A7 1 1004 0 0.000 0.00 0.00000 . 1 7 1 0.0000 0.0000 0.0000 0.0000 0.00000 2 1004 0 0.281 6.75 1.73510 0 0 7 0 2.1139 2.1139 2.1139 2.1139 1.73510 3 1004 0 8.281 198.75 0.48713 0 0 7 0 2.1139 2.1139 2.1139 2.1139 0.48713 4 1004 0 29.281 702.75 1.23790 0 0 7 0 2.1139 2.1139 2.1139 2.1139 1.23790 5 1004 0 78.283 1878.80 1.67740 0 0 7 0 2.1139 2.1139 2.1139 2.1139 1.67740 6 1004 0 176.283 4230.80 1.99880 0 0 7 0 2.1139 2.1139 2.1139 2.1139 1.99880 7 1004 0 260.283 6246.80 2.07540 0 0 7 0 2.1139 2.1139 2.1139 2.1139 2.07540 8 1004 12 364.283 8742.80 1.24820 0 0 7 0 2.1139 2.1139 2.1139 2.1139 1.24820 9 1004 12 371.283 8910.80 0.67630 0 0 7 0 2.1139 2.1139 2.1139 2.1139 0.67630 10 1004 12 378.283 9078.80 1.30550 0 0 7 0 2.1139 2.1139 2.1139 2.1139 1.30550 11 1004 12 385.283 9246.80 1.66210 0 0 7 0 2.1139 2.1139 2.1139 2.1139 1.66210 12 1004 12 393.283 9438.80 1.75560 0 0 7 0 2.1139 2.1139 2.1139 2.1139 1.75560 13 1004 12 456.292 10951.00 1.94940 0 0 7 0 2.1139 2.1139 2.1139 2.1139 1.94940 14 1004 12 561.292 13471.00 2.08330 0 0 7 0 2.1139 2.1139 2.1139 2.1139 2.08330 15 1004 12 638.292 15319.00 2.10590 0 0 7 0 2.1139 2.1139 2.1139 2.1139 2.10590 16 1004 24 728.292 17479.00 1.52450 0 0 7 0 2.1139 2.1139 2.1139 2.1139 1.52450 17 1004 24 737.292 17695.00 0.60742 0 0 7 0 2.1139 2.1139 2.1139 2.1139 0.60742 18 1004 24 743.292 17839.00 0.89626 0 0 7 0 2.1139 2.1139 2.1139 2.1139 0.89626 19 1004 24 750.292 18007.00 1.14560 0 0 7 0 2.1139 2.1139 2.1139 2.1139 1.14560 20 1004 24 757.292 18175.00 1.33590 0 0 7 0 2.1139 2.1139 2.1139 2.1139 1.33590 21 1004 24 820.292 19687.00 1.82800 0 0 7 0 2.1139 2.1139 2.1139 2.1139 1.82800 22 1006 0 0.000 0.00 0.00000 . 1 7 1 2.1139 2.1139 2.1139 2.1139 1.82800 23 1006 0 0.291 6.98 1.15610 0 0 7 0 2.1182 2.1182 2.1182 2.1182 1.15610 24 1006 0 8.291 198.98 0.33814 0 0 7 0 2.1182 2.1182 2.1182 2.1182 0.33814 25 1006 0 29.291 702.98 0.90461 0 0 7 0 2.1182 2.1182 2.1182 2.1182 0.90461 26 1006 0 80.292 1927.00 1.35550 0 0 7 0 2.1182 2.1182 2.1182 2.1182 1.35550 27 1006 0 150.292 3607.00 1.87420 0 0 7 0 2.1182 2.1182 2.1182 2.1182 1.87420 28 1006 0 246.292 5911.00 2.07530 0 0 7 0 2.1182 2.1182 2.1182 2.1182 2.07530 29 1006 12 364.292 8743.00 1.43560 0 0 7 0 2.1182 2.1182 2.1182 2.1182 1.43560 30 1006 12 371.292 8911.00 0.39132 0 0 7 0 2.1182 2.1182 2.1182 2.1182 0.39132 Now of course when I try to actually fit this data set, the model crashes due to infinite objective function. Can anyone see what I am doing wrong here? Thanks alot in advance, Pete Bonate Peter L. Bonate, PhD Genzyme Corporation Director, Pharmacokinetics 4545 Horizon Hill Blvd San Antonio, TX 78229 USA peter.bonate@genzyme.com phone: 210-949-8662 fax: 210-949-8219