Dear NONMEM Users
I’m struggling with a pharmacodynamic model for the intravenous anaesthetic,
propofol and I would really appreciate some opinions on what might be going
wrong. I have taken a sequential approach to the PK-PD modelling. PK are
described using a 3 compartment mamillary model. Bispectral Index (BIS), an
EEG derivative, was used as an effect measure. Drug was administered
intravenously (2mg/kg propofol over 1 minute) to healthy volunteers (n=6). BIS
was recorded every 15 seconds prior to drug administration and for about an
hour afterwards. BIS has a value of around 100 in an awake individual, while a
value of 40-60 indicates anaesthesia.
The data are pretty clean so I don’t understand why I’m having such difficulty.
I’ve modelled much noisier data generated with a second sedative-hypnotic drug
from this same group of patients (cross-over study) with fewer problems.
However, for this data set I have yet to produce a single run that minimises
successfully without it having a final estimate at the lower boundary for Emax
(doesn’t seem to matter how low I set the boundary). The observed Emax is
pretty low so an estimate of 20-30 wouldn’t be too unrealistic but if I set a
lower bound of -10 (NB this was just to prove the problem, I wouldn’t
ordinarily set a negative bound), NONMEM is happy to minimise with an Emax
value of -9.9. I’m using NONMEM 6 v2, FOCE, additive error. I’ve tried
additive, constant CV and log error models (with transformed data), same
problem with all.
When I model data from each subject individually, all but one (5/6) also
minimises at the lower bound for Emax so I don’t think data from anyone one
individual is causing the problem. I’ve checked for gross errors (dosing, PK
parameters). I’ve tried running with FO and the result of that is that
estimates for Emax sit on the upper boundary, rather than the lower one and the
models are strongly over-predicting.
Really hoping that I’ve not overlooked something very obvious here :-S I’ve
attached an example control stream but not the data due to limitation on the
size this e-mail could be (very happy to e-mail the data directly to anyone who
is interested in taking a look at it).
With all best wishes and very many thanks! :-)
Ann
_______________________________________________________________________
Ann Rigby-Jones PhD MRSC
Research Fellow in Pharmacokinetics & Pharmacodynamics
Peninsula College of Medicine & Dentistry
Plymouth, UK
_______________________________________________________________________
$PROB propofol PD
$INPUT ID PER TIME DV AMT RATE EVID V1 K10 K12 K21 K13 K31
$DATA BIS_Step_3_Propofol_smth.CSV IGNORE=#
$SUBROUTINES ADVAN6 TOL=3
$MODEL
COMP(CENTRAL, DEFDOSE, DEFOBS)
COMP(PERIPH1)
COMP(PERIPH2)
COMP(EFFECT)
$PK
EMAX=THETA(1)*EXP(ETA(1)) ; maximum response
E0=THETA(2)*EXP(ETA(2)) ; baseline
C50=THETA(3)*EXP(ETA(3)) ; concentration associated with 50% peak effect
GAM=THETA(4)*EXP(ETA(4)) ; gamma
K41=THETA(5)*EXP(ETA(5)) ;ke0
V4=0.00001
K14=V4*K41/V1
$DES
DADT(1)=A(2)*K21+A(3)*K31+K41*A(4)-A(1)*(K10+K12+K13+K14)
DADT(2)=A(1)*K12-A(2)*K21
DADT(3)=A(1)*K13-A(3)*K31
DADT(4)=A(1)*K14-A(4)*K41
$ERROR
CON=A(4)/V4
IF (CON.EQ.0) CON=0.0000001
TY=E0+(EMAX-E0)*(CON**GAM)/(C50**GAM+CON**GAM); SIGMOID EMAX MODEL
Y=TY + ERR(1)
W=TY
IPRED=TY
;IF(IPRED.LT.0.1) IPRED=0.1
IRES=DV-IPRED
IWRES=IRES/W
$THETA
(15,45,60) ;EMAX
(90, 98, 100 ) ;E0
(1000,2500, 8000) ;C50
(1,4, 10) ;GAMMA
(0.0001,0.2, 3) ;K41(KE0)
$SIGMA (15)
$OMEGA (0 FIX) ;EMAX
$OMEGA (0 FIX) ;E0
$OMEGA (0.01) ;C50
$OMEGA (0 FIX) ;Gamma
$OMEGA (0.01) ;KeO
$ESTIMATION METHOD=1 NOABORT MAXEVAL=9999 PRINT=5 SIGDIG=3 ;POSTHOC;INTERACTION
$COV PRINT=E
$TABLE ID TIME DV RES WRES IWRES IRES PRED IPRED EVID
ONEHEADER NOPRINT FILE=sdtab100
$TABLE ID C50 K41 EMAX E0 GAM
ONEHEADER NOPRINT FILE=patab100
PD modelling problem - Emax at lower bound?
Ann,
When you do sequential PK-PD, do not touch PK portion, it should be fixed. Thus, this is the PK:
DADT(1)=A(2)*K21+A(3)*K31-A(1)*(K10+K12+K13)
DADT(2)=A(1)*K12-A(2)*K21
DADT(3)=A(1)*K13-A(3)*K31
This is propofol concentration:
C=A(1)/V1
This is effect compartment:
DADT(4)=KE0*(C-A(4))
This is EFF model:
CONG=0
IF(C.GT.0) CONG=C**GAM
EFF=E0+(EMAX-E0)*CONG/(C50**GAM+CONG); SIGMOID EMAX MODEL
Now, the main problem, Y should be 100 minus effect, not the effect, that is why your EMAX tries to be negative. BIS vary from 0 to 100 with 100 being the baseline.
Y=100-EFF + ERR(1)
Initial conditions should be something like
(0,50,100) ;EMAX
(0, 2, 100) ;E0
Let me know of you have problem with this code, it should work.
Leonid
--------------------------------------
Leonid Gibiansky, Ph.D.
President, QuantPharm LLC
web: www.quantpharm.com
e-mail: LGibiansky at quantpharm.com
tel: (301) 767 5566
Ann Rigby-Jones wrote:
> Dear NONMEM Users
>
> I’m struggling with a pharmacodynamic model for the intravenous anaesthetic,
> propofol and I would really appreciate some opinions on what might be going
> wrong. I have taken a sequential approach to the PK-PD modelling. PK are
> described using a 3 compartment mamillary model. Bispectral Index (BIS), an
> EEG derivative, was used as an effect measure. Drug was administered
> intravenously (2mg/kg propofol over 1 minute) to healthy volunteers (n=6). BIS
> was recorded every 15 seconds prior to drug administration and for about an
> hour afterwards. BIS has a value of around 100 in an awake individual, while a
> value of 40-60 indicates anaesthesia.
>
> The data are pretty clean so I don’t understand why I’m having such difficulty.
> I’ve modelled much noisier data generated with a second sedative-hypnotic drug
> from this same group of patients (cross-over study) with fewer problems.
> However, for this data set I have yet to produce a single run that minimises
> successfully without it having a final estimate at the lower boundary for Emax
> (doesn’t seem to matter how low I set the boundary). The observed Emax is
> pretty low so an estimate of 20-30 wouldn’t be too unrealistic but if I set a
> lower bound of -10 (NB this was just to prove the problem, I wouldn’t
> ordinarily set a negative bound), NONMEM is happy to minimise with an Emax
> value of -9.9. I’m using NONMEM 6 v2, FOCE, additive error. I’ve tried
> additive, constant CV and log error models (with transformed data), same
> problem with all.
>
> When I model data from each subject individually, all but one (5/6) also
> minimises at the lower bound for Emax so I don’t think data from anyone one
> individual is causing the problem. I’ve checked for gross errors (dosing, PK
> parameters). I’ve tried running with FO and the result of that is that
> estimates for Emax sit on the upper boundary, rather than the lower one and the
> models are strongly over-predicting.
>
> Really hoping that I’ve not overlooked something very obvious here :-S I’ve attached an example control stream but not the data due to limitation on the size this e-mail could be (very happy to e-mail the data directly to anyone who is interested in taking a look at it).
>
> With all best wishes and very many thanks! :-)
>
> Ann
> _______________________________________________________________________
> Ann Rigby-Jones PhD MRSC
> Research Fellow in Pharmacokinetics & Pharmacodynamics
>
> Peninsula College of Medicine & Dentistry
> Plymouth, UK
> _______________________________________________________________________
>
> $PROB propofol PD
> $INPUT ID PER TIME DV AMT RATE EVID V1 K10 K12 K21 K13 K31
> $DATA BIS_Step_3_Propofol_smth.CSV IGNORE=#
> $SUBROUTINES ADVAN6 TOL=3
>
> $MODEL COMP(CENTRAL, DEFDOSE, DEFOBS)
>
> COMP(PERIPH1)
> COMP(PERIPH2)
> COMP(EFFECT)
>
> $PK
>
> EMAX=THETA(1)*EXP(ETA(1)) ; maximum response
>
> E0=THETA(2)*EXP(ETA(2)) ; baseline C50=THETA(3)*EXP(ETA(3)) ; concentration associated with 50% peak effect GAM=THETA(4)*EXP(ETA(4)) ; gamma K41=THETA(5)*EXP(ETA(5)) ;ke0
>
> V4=0.00001
> K14=V4*K41/V1
>
> $DES
> DADT(1)=A(2)*K21+A(3)*K31+K41*A(4)-A(1)*(K10+K12+K13+K14)
> DADT(2)=A(1)*K12-A(2)*K21
> DADT(3)=A(1)*K13-A(3)*K31
> DADT(4)=A(1)*K14-A(4)*K41
>
> $ERROR CON=A(4)/V4 IF (CON.EQ.0) CON=0.0000001
>
> TY=E0+(EMAX-E0)*(CON**GAM)/(C50**GAM+CON**GAM); SIGMOID EMAX MODEL
>
> Y=TY + ERR(1)
> W=TY
> IPRED=TY
> ;IF(IPRED.LT.0.1) IPRED=0.1
> IRES=DV-IPRED
> IWRES=IRES/W
>
> $THETA
> (15,45,60) ;EMAX
> (90, 98, 100 ) ;E0
> (1000,2500, 8000) ;C50
> (1,4, 10) ;GAMMA
> (0.0001,0.2, 3) ;K41(KE0)
>
> $SIGMA (15)
>
> $OMEGA (0 FIX) ;EMAX
> $OMEGA (0 FIX) ;E0
> $OMEGA (0.01) ;C50
> $OMEGA (0 FIX) ;Gamma
> $OMEGA (0.01) ;KeO
>
> $ESTIMATION METHOD=1 NOABORT MAXEVAL=9999 PRINT=5 SIGDIG=3 ;POSTHOC;INTERACTION
> $COV PRINT=E
> $TABLE ID TIME DV RES WRES IWRES IRES PRED IPRED EVID
> ONEHEADER NOPRINT FILE=sdtab100
>
> $TABLE ID C50 K41 EMAX E0 GAM ONEHEADER NOPRINT FILE=patab100
Ann,
Are you sure you have observed Emax in your data? The physician attending a
recently deceased American singer probably has a better insight than myself,
but it seems clear that high enough doses of propofol are fatal (and which
point presumably BIS score is 0?). Why don't you try having Emax as a
fractional decrease (constrained between 0 and 1) in E0. If Emax goes
towards 1 (i.e. a BIS of 0) then you could remove it from the model, and
just have the maximum drop in BIS to be E0. If you do estimate an Emax in
this way, I suggest a LOGIT transformation of the IIV.
Good luck,
Joe
Quoted reply history
-----Original Message-----
From: [email protected] [mailto:[email protected]] On
Behalf Of Ann Rigby-Jones
Sent: den 2 oktober 2009 15:20
To: '[email protected]'
Cc: 'Bachman, William'
Subject: [NMusers] PD modelling problem - Emax at lower bound?
Dear NONMEM Users
I'm struggling with a pharmacodynamic model for the intravenous anaesthetic,
propofol and I would really appreciate some opinions on what might be going
wrong. I have taken a sequential approach to the PK-PD modelling. PK are
described using a 3 compartment mamillary model. Bispectral Index (BIS), an
EEG derivative, was used as an effect measure. Drug was administered
intravenously (2mg/kg propofol over 1 minute) to healthy volunteers (n=6).
BIS was recorded every 15 seconds prior to drug administration and for about
an hour afterwards. BIS has a value of around 100 in an awake individual,
while a value of 40-60 indicates anaesthesia.
The data are pretty clean so I don't understand why I'm having such
difficulty. I've modelled much noisier data generated with a second
sedative-hypnotic drug from this same group of patients (cross-over study)
with fewer problems. However, for this data set I have yet to produce a
single run that minimises successfully without it having a final estimate at
the lower boundary for Emax (doesn't seem to matter how low I set the
boundary). The observed Emax is pretty low so an estimate of 20-30 wouldn't
be too unrealistic but if I set a lower bound of -10 (NB this was just to
prove the problem, I wouldn't ordinarily set a negative bound), NONMEM is
happy to minimise with an Emax value of -9.9. I'm using NONMEM 6 v2, FOCE,
additive error. I've tried additive, constant CV and log error models (with
transformed data), same problem with all.
When I model data from each subject individually, all but one (5/6) also
minimises at the lower bound for Emax so I don't think data from anyone one
individual is causing the problem. I've checked for gross errors (dosing,
PK parameters). I've tried running with FO and the result of that is that
estimates for Emax sit on the upper boundary, rather than the lower one and
the models are strongly over-predicting.
Really hoping that I've not overlooked something very obvious here :-S I've
attached an example control stream but not the data due to limitation on the
size this e-mail could be (very happy to e-mail the data directly to anyone
who is interested in taking a look at it).
With all best wishes and very many thanks! :-)
Ann
_______________________________________________________________________
Ann Rigby-Jones PhD MRSC
Research Fellow in Pharmacokinetics & Pharmacodynamics
Peninsula College of Medicine & Dentistry
Plymouth, UK
_______________________________________________________________________
$PROB propofol PD
$INPUT ID PER TIME DV AMT RATE EVID V1 K10 K12 K21 K13 K31
$DATA BIS_Step_3_Propofol_smth.CSV IGNORE=#
$SUBROUTINES ADVAN6 TOL=3
$MODEL
COMP(CENTRAL, DEFDOSE, DEFOBS)
COMP(PERIPH1)
COMP(PERIPH2)
COMP(EFFECT)
$PK
EMAX=THETA(1)*EXP(ETA(1)) ; maximum response
E0=THETA(2)*EXP(ETA(2)) ; baseline
C50=THETA(3)*EXP(ETA(3)) ; concentration associated with 50% peak effect
GAM=THETA(4)*EXP(ETA(4)) ; gamma
K41=THETA(5)*EXP(ETA(5)) ;ke0
V4=0.00001
K14=V4*K41/V1
$DES
DADT(1)=A(2)*K21+A(3)*K31+K41*A(4)-A(1)*(K10+K12+K13+K14)
DADT(2)=A(1)*K12-A(2)*K21
DADT(3)=A(1)*K13-A(3)*K31
DADT(4)=A(1)*K14-A(4)*K41
$ERROR
CON=A(4)/V4
IF (CON.EQ.0) CON=0.0000001
TY=E0+(EMAX-E0)*(CON**GAM)/(C50**GAM+CON**GAM); SIGMOID EMAX MODEL
Y=TY + ERR(1)
W=TY
IPRED=TY
;IF(IPRED.LT.0.1) IPRED=0.1
IRES=DV-IPRED
IWRES=IRES/W
$THETA
(15,45,60) ;EMAX
(90, 98, 100 ) ;E0
(1000,2500, 8000) ;C50
(1,4, 10) ;GAMMA
(0.0001,0.2, 3) ;K41(KE0)
$SIGMA (15)
$OMEGA (0 FIX) ;EMAX
$OMEGA (0 FIX) ;E0
$OMEGA (0.01) ;C50
$OMEGA (0 FIX) ;Gamma
$OMEGA (0.01) ;KeO
$ESTIMATION METHOD=1 NOABORT MAXEVAL=9999 PRINT=5 SIGDIG=3
;POSTHOC;INTERACTION
$COV PRINT=E
$TABLE ID TIME DV RES WRES IWRES IRES PRED IPRED EVID
ONEHEADER NOPRINT FILE=sdtab100
$TABLE ID C50 K41 EMAX E0 GAM
ONEHEADER NOPRINT FILE=patab100
Dear Ann,
I would recommend you to investigate the possibility to model the BIS and
propofol data using a two-effect site model instead of just a single effect
site model. We presented a two-effect site model for BIS and propofol at PAGE
this year (Reference: PAGE 18 (2009) Abstr 1590
[www.page-meeting.org/?abstract=1590]). This approach gave substantially better
predictions than just using a single effect site model.
Best regards,
Ulrika Simonsson
Assoc Prof
The Pharmacometrics Research Group
Uppsala University, Sweden
Quoted reply history
-----Original Message-----
From: [email protected] [mailto:[email protected]] On
Behalf Of Leonid Gibiansky
Sent: den 2 oktober 2009 16:11
To: Ann Rigby-Jones
Cc: '[email protected]'
Subject: Re: [NMusers] PD modelling problem - Emax at lower bound?
Ann,
When you do sequential PK-PD, do not touch PK portion, it should be
fixed. Thus, this is the PK:
DADT(1)=A(2)*K21+A(3)*K31-A(1)*(K10+K12+K13)
DADT(2)=A(1)*K12-A(2)*K21
DADT(3)=A(1)*K13-A(3)*K31
This is propofol concentration:
C=A(1)/V1
This is effect compartment:
DADT(4)=KE0*(C-A(4))
This is EFF model:
CONG=0
IF(C.GT.0) CONG=C**GAM
EFF=E0+(EMAX-E0)*CONG/(C50**GAM+CONG); SIGMOID EMAX MODEL
Now, the main problem, Y should be 100 minus effect, not the effect,
that is why your EMAX tries to be negative. BIS vary from 0 to 100 with
100 being the baseline.
Y=100-EFF + ERR(1)
Initial conditions should be something like
(0,50,100) ;EMAX
(0, 2, 100) ;E0
Let me know of you have problem with this code, it should work.
Leonid
--------------------------------------
Leonid Gibiansky, Ph.D.
President, QuantPharm LLC
web: www.quantpharm.com
e-mail: LGibiansky at quantpharm.com
tel: (301) 767 5566
Ann Rigby-Jones wrote:
> Dear NONMEM Users
>
> I’m struggling with a pharmacodynamic model for the intravenous anaesthetic,
> propofol and I would really appreciate some opinions on what might be going
> wrong. I have taken a sequential approach to the PK-PD modelling. PK are
> described using a 3 compartment mamillary model. Bispectral Index (BIS), an
> EEG derivative, was used as an effect measure. Drug was administered
> intravenously (2mg/kg propofol over 1 minute) to healthy volunteers (n=6).
> BIS was recorded every 15 seconds prior to drug administration and for about
> an hour afterwards. BIS has a value of around 100 in an awake individual,
> while a value of 40-60 indicates anaesthesia.
>
> The data are pretty clean so I don’t understand why I’m having such
> difficulty. I’ve modelled much noisier data generated with a second
> sedative-hypnotic drug from this same group of patients (cross-over study)
> with fewer problems. However, for this data set I have yet to produce a
> single run that minimises successfully without it having a final estimate at
> the lower boundary for Emax (doesn’t seem to matter how low I set the
> boundary). The observed Emax is pretty low so an estimate of 20-30 wouldn’t
> be too unrealistic but if I set a lower bound of -10 (NB this was just to
> prove the problem, I wouldn’t ordinarily set a negative bound), NONMEM is
> happy to minimise with an Emax value of -9.9. I’m using NONMEM 6 v2, FOCE,
> additive error. I’ve tried additive, constant CV and log error models (with
> transformed data), same problem with all.
>
> When I model data from each subject individually, all but one (5/6) also
> minimises at the lower bound for Emax so I don’t think data from anyone one
> individual is causing the problem. I’ve checked for gross errors (dosing, PK
> parameters). I’ve tried running with FO and the result of that is that
> estimates for Emax sit on the upper boundary, rather than the lower one and
> the models are strongly over-predicting.
>
> Really hoping that I’ve not overlooked something very obvious here :-S I’ve
> attached an example control stream but not the data due to limitation on the
> size this e-mail could be (very happy to e-mail the data directly to anyone
> who is interested in taking a look at it).
>
> With all best wishes and very many thanks! :-)
>
> Ann
> _______________________________________________________________________
> Ann Rigby-Jones PhD MRSC
> Research Fellow in Pharmacokinetics & Pharmacodynamics
>
> Peninsula College of Medicine & Dentistry
> Plymouth, UK
> _______________________________________________________________________
>
> $PROB propofol PD
> $INPUT ID PER TIME DV AMT RATE EVID V1 K10 K12 K21 K13 K31
> $DATA BIS_Step_3_Propofol_smth.CSV IGNORE=#
> $SUBROUTINES ADVAN6 TOL=3
> $MODEL
> COMP(CENTRAL, DEFDOSE, DEFOBS)
> COMP(PERIPH1)
> COMP(PERIPH2)
> COMP(EFFECT)
>
> $PK
>
> EMAX=THETA(1)*EXP(ETA(1)) ; maximum response
> E0=THETA(2)*EXP(ETA(2)) ; baseline
> C50=THETA(3)*EXP(ETA(3)) ; concentration associated with 50% peak effect
> GAM=THETA(4)*EXP(ETA(4)) ; gamma
> K41=THETA(5)*EXP(ETA(5)) ;ke0
>
>
> V4=0.00001
> K14=V4*K41/V1
>
> $DES
> DADT(1)=A(2)*K21+A(3)*K31+K41*A(4)-A(1)*(K10+K12+K13+K14)
> DADT(2)=A(1)*K12-A(2)*K21
> DADT(3)=A(1)*K13-A(3)*K31
> DADT(4)=A(1)*K14-A(4)*K41
>
> $ERROR
> CON=A(4)/V4
> IF (CON.EQ.0) CON=0.0000001
>
> TY=E0+(EMAX-E0)*(CON**GAM)/(C50**GAM+CON**GAM); SIGMOID EMAX MODEL
>
> Y=TY + ERR(1)
> W=TY
> IPRED=TY
> ;IF(IPRED.LT.0.1) IPRED=0.1
> IRES=DV-IPRED
> IWRES=IRES/W
>
>
> $THETA
> (15,45,60) ;EMAX
> (90, 98, 100 ) ;E0
> (1000,2500, 8000) ;C50
> (1,4, 10) ;GAMMA
> (0.0001,0.2, 3) ;K41(KE0)
>
> $SIGMA (15)
>
> $OMEGA (0 FIX) ;EMAX
> $OMEGA (0 FIX) ;E0
> $OMEGA (0.01) ;C50
> $OMEGA (0 FIX) ;Gamma
> $OMEGA (0.01) ;KeO
>
> $ESTIMATION METHOD=1 NOABORT MAXEVAL=9999 PRINT=5 SIGDIG=3
> ;POSTHOC;INTERACTION
> $COV PRINT=E
> $TABLE ID TIME DV RES WRES IWRES IRES PRED IPRED EVID
> ONEHEADER NOPRINT FILE=sdtab100
> $TABLE ID C50 K41 EMAX E0 GAM
> ONEHEADER NOPRINT FILE=patab100
>