From: "Eleveld, DJ"
Subject: [NMusers] Modelling twitch potentiation in NONMEM, problems with $DES 
Date: Thu, March 10, 2005 6:22 am 

Hello everyone, 

I am wondering if anyone can give me advice or insight into using
NONMEM for modelling twitch potentiation along with a pharmacodynamic
model.  I am trying to use NONMEM to model the method described in "Model
to describe the degree of twitch potentiation during neuromuscular
monitoring, D. J. Eleveld, A. F. Kopman, J. H. Proost, J. M. K. H. Wierda,
British journal of Anaesthesia, Br J Anaesth 2004; 92: 373-80".

The basic idea is that the degree of twitch potentiation increases as a
result of a twitch and decreases exponentially in the inter-twitch period.
I have been using a compartment to model the degree of potentiation.  Because
I need to increase the degree of potentiation step-wise immediately after
each twitch observation, I have used another compartment to act as a
counter with the twitch observation count.  

The fitting is related to the standard LN-PD model fitting for muscle
relaxants.  It is slightly more complex because the twitch height calculation
has an additional factor varies dependant on an 'internal' variable, the NMB.
What I mean by 'internal' here is that NMB is dependent on Keo, which I am
also fitting.  It seems to be logical to model the degree of potentiation
with a compartment, but then I have to make sure this compartment gets a
'dose' at every twitch observation.  I can't calculate the 'dose' size
before hand because it depends on the NMB, which depends on the Keo, which
I am also fitting.  So I think I have to incorporate the 'dose' into the
potentiation compartment within the $DES code.  This doesn't seem to be the
usual way of writing $DES code, but I can't think of any other way to get it to work.

There is the slight added complication that there are some obviously incorrect
twitches removed from the data file, but these twitches cannot simply be
ignored as they do include a twitch and therefore influence the time course
of potentiation.  These have been marked as 999 in the data file, a number
that otherwise does not occur.

My control file looks like: 

$PROB  Potentiation LN-PD fitting 
$DATA  data2.prn 
$MODEL COMP(EFFECT)    ; NMB effect compartment 
       COMP(POTENT)    ; Potentiation compartment 
       COMP(HACK)      ; Ugly (but necessary) twitch counter 
    KEO=THETA(1)*EXP(ETA(1))   ; Link Keo 
    EC50=THETA(2)*EXP(ETA(2))  ; PD Ec50 
    GAMM=THETA(3)*EXP(ETA(3))  ; PD Gamma 
    SCAL=THETA(4)+ETA(4)       ; Potentiation model Kscale 
    POTR=ABS(THETA(5)+ETA(5))  ; Potentiation model r must be +ve 
    POTK=THETA(6)*EXP(ETA(6))  ; Potentiation model k 
    DADT(1)=KEO*(PCON-A(1))    ; Effect site concentration decay 
    PDT1=A(1)**GAMM            ; Calculate NMB 
    PINC=POTR*(1-NMB)          ; Increase potentiation 
    IF (A(3).EQ.PSTP) PINC=0   ; only once per twitch 
    DADT(2)=-POTK*A(2)+PINC    ; The degree of potentiation (P) 
    DADT(3)=-A(3)+PSTP         ; Keep twicth count flag correct 
    EDT1=A(1)**GAMM            ; NMB must be diff names than $DES, sigh.... 
    Y=SCAL*(BASE/100)*(1-ENMB)+ERR(1) ; Predicted twitch response 
    IF (DV.EQ.999) Y=999+ERR(1); Missing twitches are fit almost excatly 
$THETA (0.05,,0.25)(1000.,,2500.)(2.,,6.)(90.,,120.) 
$OMEGA 0.5 0.5 0.5 10 0.01 0.5 
$SIGMA 20 

This way to model increases in potentiation seems logical, if a bit awkward.  The amount
in compartment 2 is the degree of twitch potentiation and compartment 3 functions as a
'counter' along with the PSTP record (different for each observation) to make sure
potentiation only increases once per observation even if the $DES code is called more
than once per observation.  

Possibly my assumptions about how the $DES code is called are incorrect.  Does anyone
see any obvious problems with my control file or reasons why it might not work as expected?

Thank you, 

Doug Eleveld