NONMEM versus Phoenix NLME: Discontinuity in differential equation

Phoenix Certara Forums NLME-NONMEM model comparisons
1

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:P1exp(-khd1(tsld-0.01)))

HD2=(tsld<1?0:P2exp(-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

2

you can do conditional ODE and compute time since last intervention refer to :
https://support.certara.com/forums/topic/811-ifnewindne2-equivalent-for-average-concentration-calcuation-on-the-fly/?hl=tprev#entry3326

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

3

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

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:P1exp(-khd1(tsld-0.01)))
HD2=(tsld<1?0:P2exp(-khd2(tsld-1)))

error(QTcEps = 2.49018)
observe(QTcObs = QTc + QTcEps)

#Drug effect model

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

kin=QTc0kout
Edrug=slopeC

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

#structural parameters

stparm(QTc0 = tvQTc0 * exp(nQTc0))
stparm(Amp1= tvAmp1exp(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))
Example.phxproj (324 KB)

5

Example.phxproj (324 KB)