Jump to content


Photo

NONMEM versus Phoenix NLME: Discontinuity in differential equation


  • Please log in to reply
4 replies to this topic

#1 Pradeepsharma

Pradeepsharma

    Newbie

  • Members
  • Pip
  • 9 posts

Posted 12 February 2020 - 12:30 AM

Hi,

 

I was trying to reproduce model from reference Snelder et al (2014) British Journal of Pharmacology 171 5076–5092.

 

It is indirect response model for heart rate changes in dog including (i) feedback mechanisms and (ii) alterations in heart rate changes due to manual handling.

 

Equations are :

 

dHR/dt =Kin (1+ CR).(1-FB.MAP).(1+HD)-Kout.HR

 

where

HR=heart rate

CR is circadian rhythm = amp.cos(2p (t+hor)/24)

FB= feedback component

MAP= maximum arterial pressure

Kin and kout are usual constants

HD= handling effect=P. exp[-kHD.(t-thd) when t>thd

 

Authors had used NONMEM to work out this model.

 

When I tried to use it in Phoenix NLME, it failed to fit the data and I think it is due to 'discontinuity' introduced by 'time event t-thd' in handling effect. 

 

Can anyone suggest how to overcome this using sequence statement or alternative PML syntax.  I want to use following 2 handling effects to model QTc change in indirect resposne model:

 

 

         Handling effects: (repeat of handling effect every 24 hour modelled using 'tsld')

         

          deriv(tsld=1)

          dosepoint(Aa,doafter={tsld=0})

 

          HD1=(tsld<0.01?0:P1*exp(-khd1*(tsld-0.01)))

          HD2=(tsld<1?0:P2*exp(-khd2*(tsld-1)))

 

        Circadian rhythm:

        CR = Amp1 * cos (2 * 3.1416 * (t - Phi1)/24)       

 

          Indirect response model for QTc change

           sequence{QTc0=(kin)/(kout)}

        

            Edrug=slope*C

                          

            deriv(QTc = ((kin)*(1+CR)*(1=HD1)(1+HD2+Edrug))-kout*QTc)

 

 

Best wishes

Pradeep

 

 


  • AlexKawn and RonaldAxon like this

#2 smouksassi1

smouksassi1

    Advanced Member

  • Val_Members
  • PipPipPip
  • 232 posts
  • LocationMontreal

Posted 12 February 2020 - 01:47 PM

you can do conditional ODE and compute time since last intervention refer to :
https://support.cert...tprev#entry3326

 

unless you provide a project and clean code it will be not possible to help



#3 smouksassi1

smouksassi1

    Advanced Member

  • Val_Members
  • PipPipPip
  • 232 posts
  • LocationMontreal

Posted 12 February 2020 - 01:59 PM

dosepoint(A1, idosevar = A1Dose, infdosevar = A1InfDose, infratevar = A1InfRate)
deriv(A1 = t > 5 ? - kHD*A1  : 0)
Agraph = A1
dosepoint(A1, idosevar = A1Dose, infdosevar = A1InfDose, infratevar = A1InfRate)
error(CEps = 1)
observe(CObs = A1 + CEps)
stparm(kHD = tvkHD * exp(nkHD))
fixef(tvkHD = c(, 1, ))
ranef(diag(nkHD) = c(1))
}


#4 Pradeepsharma

Pradeepsharma

    Newbie

  • Members
  • Pip
  • 9 posts

Posted 13 February 2020 - 01:18 AM

Many thanks for advise. Please see attached the project file (Example.phxproj) with detailed model.

 

My query is: - Handling effects (HD1 and HD2) have a discrete events and bring discontinuity in main drug effect differential equation. How to overcome this with alternative PML syntax?

 

The model PML code is essentially as below.

 

test(){
covariate(nV,nCl,nKa,nV2)
deriv(tsld=1)
    dosepoint(Aa,doafter={tsld=0})
 
#PK Model
deriv(Aa = - Ka * Aa)
deriv(A1 = Ka * Aa - Cl * C - Cl2 * (C - C2))
deriv(A2 = Cl2 * (C - C2))
dosepoint(Aa)
C = A1 / V
C2 = A2 / V2
error(CEps(freeze) = 0.374911)
observe(CObs = C * (1 + CEps))
 
#Baseline model
CR = Amp1 * cos (2 * 3.1416 * (t - Phi1)/24)     
       
        HD1=(tsld<0.01?0:P1*exp(-khd1*(tsld-0.01)))
            HD2=(tsld<1?0:P2*exp(-khd2*(tsld-1)))
                                       
        error(QTcEps = 2.49018)
            observe(QTcObs = QTc + QTcEps)
            
       
#Drug effect model
 
sequence{QTc0=(kin)/(kout)}
       
kin=QTc0*kout
    Edrug=slope*C
                         
    deriv(QTc = ((kin)*(1+CR)*(1+HD1)*(1+HD2+Edrug))-kout*QTc)
 
 
#structural parameters
 
   stparm(QTc0 = tvQTc0 * exp(nQTc0)) 
   stparm(Amp1= tvAmp1*exp(nAmp1))
   stparm(Phi1= tvPhi1* exp(nPhi1))
   stparm(khd1= tvkhd1* exp(nkhd1))
   stparm(khd2= tvkhd2* exp(nkhd2))  
   stparm(kin = tvkin * exp(nkin))
   stparm(kout = tvkout * exp(nkout))         
    stparm(P1 = tvP1 * exp(nP1))
stparm(P2 = tvP2 * exp(nP2))   
stparm(slope= tvslope*exp(nslope)) 
 
 
   fixef(tvQTc0 = c(,216,))
   fixef(tvAmp1 = c(,0.1,))
   fixef(tvPhi1 = c(,1,))
   fixef(tvAmp2 = c(,0.1,))
   fixef(tvPhi2 = c(,1,))
   fixef(tvP1 = c(,-9.37923,))
   fixef(tvP2 = c(,0.131449,))
   fixef(tvkhd1 (freeze)= c(,1.69524,))
   fixef(tvkhd2 (freeze)= c(,0.0121554,))
   fixef(tvslope = c(,0.0045,))
fixef(tvkin = c(,1,))
fixef(tvkout = c(,0.051,))
 
stparm(Ka = tvKa * exp(nKa))
stparm(V = tvV * exp(nV))
stparm(V2 = tvV2 * exp(nV2))
stparm(Cl = tvCl * exp(nCl))
stparm(Cl2 = tvCl2)
fixef(tvKa (freeze)= c(, 2.13259, ))
fixef(tvV (freeze)= c(, 3.82764, ))
fixef(tvV2 (freeze)= c(, 1.68737, ))
fixef(tvCl (freeze)= c(, 0.479576, ))
fixef(tvCl2 (freeze)= c(, 0.202842, ))
 
 
ranef(diag(nslope,nke0,nkin,nkout,nPhi1,nAmp1,nQTc0,nkhd1,nkhd2,nP1,nP2) = c(1,1,1,1,1,1,1,1,1,1,1))

Edited by Pradeepsharma, 13 February 2020 - 01:32 AM.


#5 Pradeepsharma

Pradeepsharma

    Newbie

  • Members
  • Pip
  • 9 posts

Posted 13 February 2020 - 01:22 AM

Attached File  Example.phxproj   324.06KB   5 downloads Looks my previous post did not include project file. So trying another time and hopefully works.


Edited by Pradeepsharma, 13 February 2020 - 01:34 AM.





0 user(s) are reading this topic

0 members, 0 guests, 0 anonymous users