Reproduce simulation from NONMEM. Error in configsaem: covariate(s) not found: A1, A2, q, centr, Vc

[AprilCCC]

Hi nlmixr team,

I just tried to reproduce a NONMEM popPK dose simulation from observed real data with nlmixr, and when I tried to simulate it, occurred 2 errors, and under the same dose, however, the NCA result from nlmixr model-simulated data was far from the original observed data NCA result. Can anyone help me to check what is wrong with the nlmixr model? The observed real data is in the attachment. Many thanks!
observed data.xlsx

Error in .nlmixrUIModel(fun, ini, bigmodel) :
There must be at least one prediction in the model({}) block. Use ~ for predictions
Error in configsaem(model = model, data = dat, inits = inits, mcmc = .mcmc, :
covariate(s) not found: A1, A2, q, centr, Vc

The NONMEN model code showed below:
$SUBROUTINE ADVAN6 TOL=3
$MODEL NCOMP=2 COMP=(CENTRAL,DEFOBS,DEFDOSE) COMP=(PER)
$PK

TVCLS = THETA(1)(WT/70) ** THETA(7)
CLS = TVCLS EXP(ETA(1))

TVVC = THETA(2)*(WT/70) ** THETA(8)
VC = TVVC * EXP(ETA(2))

TVCLP = THETA(3)(WT/70) ** THETA(7)
CLP = TVCLP EXP(ETA(3))

TVVP = THETA(4)*(WT/70) ** THETA(8)
VP = TVVP * EXP(ETA(4))

TVVMAX = THETA(5)
VMAX = TVVMAX* EXP(ETA(5))

TVKM = THETA(6)
KM = TVKM* EXP(ETA(6))

;S2=V2

$DES

C1 = A(1)/VC
C2 = A(2)/VP

DADT(1) = -VMAXC1/(KM+C1)-(CLSC1)-(CLPC1)+(CLPC2)
DADT(2) = (CLPC1)-(CLPC2)

$ERROR

R=A(1)/VC
IPRED = R
IRES = DV-IPRED
IWRES = IRES/DV
Y = R*(1+EPS(1))+EPS(2)

$THETA 0.0223 FIX ; 1CLS,L/h
3.26 FIX ; 2VC,L
0.0433 FIX ; 3CLP,L/h
2.72 FIX ; 4VP,L
295.4 FIX ; 5VMAX,ug/h
115 FIX ; 6KM,ng/ml
0.607 FIX ; 7WT on CL
0.371 FIX ; 8WT on V
;0 FIX ;WT on VMAX
$OMEGA 0.093 FIX ; ETA_CLS
0 FIX ; ETA_VC
0 FIX ; ETA_CLP
0.87 FIX ; ETA_VP
0.161 FIX ; ETA_VAMX
0 FIX ; ETA_KM
$SIGMA 0.284 FIX ; EPS(1)
14.4 FIX ; EPS(2) ng/ml
$SIMULATION (889216570) ONLYSIM SUBPROBLEM=10

The nlmixr model code I created showed below:
cmt2sad5 <- function(){
ini({
lv <- 3.26
lcl <- 0.0223
lQ <- 0.0433 # Q inter-compartmental CL (Q)
lvp <- 2.72 # Vp peripheral volume (Vp)
lVM <- 295.4
lKM <- 115
WTcl <- 0.607
WTv <-0.371
wt.est <- 0.0
eta.v ~ 0
eta.cl ~ 0.093
eta.Q ~ 0
eta.vp ~ 0.87
eta.VM ~ 0.161
eta.KM ~ 0
EPS1 ~ 0.284
EPS2 ~ 14.1
})
model({
cl <- exp(lcl + eta.cl + wt.est * (WTcl/70))
vc <- exp(lv + eta.v + wt.est * (WTv/70))
Q <- exp(lQ + eta.Q + wt.est * (WTcl/70))
vp <- exp(lvp + eta.vp + wt.est *(WTv/70))
VM <- exp(lVM + eta.VM)
KM <- exp(lKM + eta.KM)
c1 = A1 / vc
c2 = A2 / vp
d/dt(center) = -VM * c1/(KM+c1) - (cl c1) -(q * c1) + (q * c2)
d/dt(periph) = (q * c1) + (q * c2)
##linCmt() ~ lnorm(logn.sd)
R=A1 / vc
IPRED = R
IRES = DV-IPRED
IWRES = IRES/DV
Y = R(1+EPS1)+EPS2
})
}
cmt2msad5 <- nlmixr(cmt2sad5)
##print(cmt2msad5)
cmt2fitsad5.logn <- nlmixr(cmt2msad5, sad5, "saem",
control=list(print=0),
table=tableControl(cwres=TRUE, npde=TRUE))
plot(augPred(cmt2fitsad5.logn))

##simulation with same dose, 5mg/kg in human with 61.5kg weight in avarage
ev <- eventTable(amount.units="mg", time.units="h") %>%
et(amt=5 * 61.5 * 1000) %>%
et(seq(0, 504, by=1)) %>%
et(id=1:60)

sim1 <- simulate(cmt2fitsad5.logn, events=ev)

[mattfidler]

With nlmixr you cannot simulate from the un-fitted model. The upcoming nlmixr2 allows you to simulate.

However, regardless of the method you use, you need to specify residual error with the ~ syntax. This would be something like R ~ prop(e1.sd)+add(e2.sd). You also do not have to specify IPRED, IRES, IWRES in the nlmixr estimation. They are added automatically by nlmixr/nlmixr2

[AprilCCC]

Hi Matthew, thanks for the suggesion.
I still have several questions regarding reproduce model from NONMEM code, Can you help me to check with this? Really appreciate it.

1. NONMEM code: TVCLS = THETA(1)*(WT/70)*THETA(7) CLS = TVCLS EXP(ETA(1)) $THETA 0.607 FIX ; 7WT on CL
   Is it correct way to write in the nlmirx? cl <- exp(lcl + eta.cl + log(WT/70) * WTcl)
2. NONMEM code:

    $DES
   C1 = A(1)/VC
   C2 = A(2)/VP
   DADT(1) = -VMAXC1/(KM+C1)-(CLSC1)-(CLPC1)+(CLPC2)
   DADT(2) = (CLPC1)-(CLPC2)

Is it correct way to write in the nlmixr as below?

   c1 = center / vc
   c2 = periph / vp
   d/dt(center) = -VM * c1/(KM+c1) - (cl *c1) -(Q * c1) + (Q * c2)
   d/dt(periph) = (Q * c1) - (Q * c2)

4. How to specify $SIGMA 0.284 FIX ; EPS(1). 14.4 FIX ; EPS(2) ng/ml. in the nlmixr estimation?

[mattfidler]

omega values cannot be non-positive definite. NONMEM allows 0 estimates of ETAs to remove them from the calculation. nlmixr does not currently do this for users Instead you can drop these values by piping, ie %>% vc <- exp(lv + wt.est * (WTv/70))

[AprilCCC]

Hi Matthew,

These are the PK parameters simulated from NONMEM, and I have been trying to reproduce them in nlmixr. The NONMEM model codes were shown above, but I failed to simulate by Fitting the data with nlmixr, maybe the model in the nlmixr I created was not fit with observed data.

And then I tried Population Simulations with RxODE, codes shown below, and the result looks close to the PK parameters calculated from NONMEN simulated data.

• Create RxODE model for two compartment IV infusion

mod2cmt5 <- RxODE({
  cl <- tcl * exp(eta.cl * log (WT/70) * WTcl)
  vc <- tv+log(WT/70)*WTv
  Q <- tQ * exp( log (WT/70) * WTcl)
  vp <- tvp * exp(eta.vp *  log(WT/70) * WTv)
  VM <- exp(tVM + eta.VM )
  KM <- tKM
  c1 = center / vc
  c2 = periph / vp
  d/dt(center) = -VM * c1/(KM+c1) - (cl *c1) -(Q * c1) + (Q * c2)
  d/dt(periph) = (Q * c1) - (Q * c2)
})

nsub = 12. ##simulate subjects number 
theta <- cbind(tcl=0.0223,tv=3.26,tvp=2.72,tQ=0.0433,tVM=295.4,tKM=115,
           WTcl=0.607, WTv=0.371, WT = rnorm(nsub, 63, 2))
omega <- lotri(eta.cl ~ 0.093, eta.vp ~ 0.87, eta.VM ~ 0.161)
sigma <- lotri(add.err ~ 14.4, prop.err ~ 0.284)

simulation with new dose sad5

sim = list()
for (i in 1:nsub) {
  theta1 = theta[i,]
  evsad5 <- eventTable(amount.units="mg", time.units="h")
  evsad5$add.dosing(dose=4.4 * theta1["WT"] * 1000, dur=1.1)
  evsad5$add.sampling(0:504)
  sim[[i]] = cbind(id=i,rxSolve(mod2cmt5,evsad5, theta1, omega=omega,sigma=sigma,addDosing=TRUE))
}

simsad5 = do.call("rbind",sim)
simsad5data <- as.data.frame(simsad5)
print(simsad5data)
 
sad5_conc <- PKNCAconc(data=simsad5data, formula = c1 ~ time | id )
sad5_dose <- PKNCAdose(data=simsad5data, formula = amt ~ time | id, route = "intravascular", dose.time = 0,duration = 1.1)
pk_sad5sim <- PKNCAdata(sad5_conc, sad5_dose, intervals=data.frame(start=0, end=504,
        tmax=TRUE,cmax=TRUE, thalf.eff.iv.last =TRUE, aucinf.obs=TRUE,lambda.z=TRUE,cl.obs=TRUE,vz.obs=TRUE))
result_sad5sim <- pk.nca(pk_sad5sim)
PKNCA.set.summary(name=c("thalf.eff.iv.last "),point=business.median,spread=business.range,description="median and range")
summary(result_sad5sim)

The c1 concentration NCA result shown below:

sstart end   N         cmax                      tmax                               thalf.eff.iv.last            lambda.z                aucinf.obs
     0   504    12   82000 [2.10]    2.00 [2.00, 2.00]           67.2 [65.5, 67.7]          0.0177 [0.731]       6.85e6 [2.87]
      cl.obs                     vz.obs
 0.0403 [0.770]      2.27 [1.50]

Is it correct that include simulated c1 (central concentration from two compartment model) for PKNCA calculation? The simulated Cmax and AUC data look a bit close to NONMEN simulated data, but how can I adjust the code to look closer, especially for Cl, Vz, and C336h?
Should I just exclude parameters in the RxODE model when 0 estimates of ETAs in NONMEN?

Can you help me to check the codes? Thank you.

[mattfidler]

I'm unsure. If the code is reproducing the NONMEM code, it is likely OK.

[AprilCCC]

Hi @mattfidler

For the model I created above, if I want to add more intersubject variations on c1 and cl, how should I adjust sigma? should that be close to 0 or a bigger number? Thanks.

[mattfidler]

The inter-subject variability can be changed by adding larger omega variations. The sigma adds larger residual variations