Skip to contents

Simulation with Modifications to Parameters and Model Code

Philip Delff

Simone Cassani

June 12, 2025 using NMsim 0.2.3.906

Source: vignettes/NMsim-modify-model.Rmd
NMsim-modify-model.Rmd

Introduction

This vignette introduce NMsim features that allows users to specify modifications to a NONMEM model. The NONMEM control stream can be edited directly from within the NMsim function call in the R script via the flexible interface provided by the argument modify.model. The edits can be applied to any section of the control stream, and most commonly they will affect the $PK or $PRED sections, or even the $DES or $MODEL sections when needed. Note that, since the data related sections ($INPUT,$DATA,$TABLE) are automatically handled by NMsim, it will be rarely necessary that they are edited via the modify.model argument. Additionally the NMsim function has specific arguments regarding the $TABLE sections, namely table.vars and table.options.

Notice, until NMsim 0.1.0, the modify.model argument was called list.sections, and it was characterized by a different formatting and level of flexibility. The helper functions add and overwrite were introduced with NMsim 0.1.3, to add and overwrite control stream lines, respectively. Until then, custom modification functions had to be provided like described in the section “Edit control stream using customized functions’.

This vignette illustrates some examples and applications of the use of the modify.model argument.

An initial configuration, described in the Introductory NMsim vignette add_link, might be necessary for the correct functioning of the NMsim library.

Objectives

The examples described below use the modify.model argument to modify the NONMEM control stream and to address the following pharmacometric questions:

  • Modify KA: How is the concentration-time profile affected if switching formulation (reducing absorption rate) by a range of fold-values?
  • Modify F1 and CL: What is the expected concentration-time profile in patients with a certain Drug-Drug Interaction effect on clearance and bioavailability?
  • Add ALAG: How will a dose delay of different amounts of time affect the predicted exposure?
  • Add AUC computation in control stream: How can we use NONMEM to integrate AUC during model simulation leveraging the benefits of using a sparse grid via NMsim?

Example: Change in formulation

The effect on the concentration-time profile of a change in formulation that reduces absorption rate KA is explored with the use of a scaling factor KASCALE=c(1,4,10) included to the NONMEM control stream via the modify.model argument. The absorption rate reduction is provided through the NMsim simulation data dat.sim.varka containing dosing events, simulation time steps and the value of KASCALE for each simulated patient.

First KASCALE is added to the simulation data

file.mod <- system.file("examples/nonmem/xgxr021.mod",
                        package="NMsim")
data.sim <- NMreadCsv(system.file("examples/derived/dat_sim1.csv",
                                 package="NMsim"))

#data.sim=setDT(data.sim)
# add KASCALE and copy patient info for each value of KASCALE
dat.sim.varka <- data.sim[,data.table(KASCALE=c(1,4,10)),by=data.sim] 
dat.sim.varka[,ID:=.GRP,by=.(KASCALE,ID)] # update patient IDs
setorder(dat.sim.varka,ID,TIME,EVID) # order rows

Then the NMsim function is used to run the simulation in which the modify.model argument is used to simulate a modified model with different absorption rates, scaled by the parameter KASCALE. Two different approaches are demonstrated:

  1. the effect of the scaling factor KASCALE is added at the end of the PK section in the NONMEM control stream.
simres.varka <- NMsim(file.mod=file.mod # NONMEM control stream
          ,data=dat.sim.varka # simulation data file
          ,name.sim="varka"
          ,modify.model=list(PK=add("KA=KA/KASCALE")))
  1. the specific line that defines TVKA is modified to include the effect of the scaling factor KASCALE 
simres.varka2 <- NMsim(file.mod=file.mod
          ,data=dat.sim.varka
          ,name.sim="varka2"
          ,modify.model=list(PK=overwrite("THETA(1)","THETA(1)/KASCALE"))
          )

The code below returns the plot

simres.varka=as.data.table(simres.varka)
simres.varka2=as.data.table(simres.varka2)
simres.both <- rbind(simres.varka[,method:="a. add()"],
                     simres.varka2[,method:="b. custom/sub()"]
                     )

ggplot(simres.both[EVID==2],aes(TIME,PRED,colour=factor(KASCALE)))+
  geom_line()+
  labs(colour="Fold absorption prolongation, KASCALE")+
  scale_x_continuous(breaks=seq(0,168,by=24))+
  facet_wrap(~method)+
  scale_color_manual(values=c("orange", "blue", "darkgreen"))
Concentration (`PRED`) profile as a function of time computed by `NMsim` modified models **a** (left) and **b** (right). The equivalence and robustness of the two modified models is supported by the matching results, corresponding to reduced `PRED` values for higher values of `KASCALE` (lower absorption rate).

Concentration (PRED) profile as a function of time computed by NMsim modified models a (left) and b (right). The equivalence and robustness of the two modified models is supported by the matching results, corresponding to reduced PRED values for higher values of KASCALE (lower absorption rate).

NOTE: the NMsim overwrite function prior to version 0.1.6 has a bug that was fixed in the0.1.5.902 github release; it can be installed with

library(remotes)
install_github("NMautoverse/NMsim")

Drug-Drug Interaction (DDI)

The effect of DDI on clearance (CL) and bioavailability (F1) is simulated for the following scenarios
scenario CLSCALE FSCALE CLSCALE/FSCALE
noDDI 1 1 1
DDI.1 0.5 1.2 0.42
DDI.2 0.33 1.1 0.3

First the CL and F1 data is added to the patient data

## 'Outer join'
dat.sim.DDI=data.sim[,data.table(CLSCALE=c(1,1/2,1/3)
                                ,FSCALE=c(1,1.2,1.1)
                                ,lab=c("noDDI","DDI.1","DDI.2"))
                    ,by=data.sim]
dat.sim.DDI[,ID:=.GRP,by=.(lab,ID)]
setorder(dat.sim.DDI,ID,TIME,EVID)

Then, the DDI driven change in parameters is added at the end of the PK section of the NONMEM control stream via the modify.model argument in the NMsim() function.

simres.DDI <- NMsim(file.mod=file.mod
            ,data=dat.sim.DDI
            ,name.sim="DDI"
            ,modify.model=list(PK=add("CL=CL*CLSCALE"
                                      ,"F1=FSCALE")))

The effect on the concentration-time profile is shown in the figure below

Concentration (`PRED`) profile as a function of time computed by `NMsim` modified model for different DDIs. The modified model correctly simulates (i) a higher value of `Cmax` on day 1 for higher biovalability; (ii) a higher `PRED` value at steady state for lower apparent clearance effect `CLSCALE/FSCALE` values.

Concentration (PRED) profile as a function of time computed by NMsim modified model for different DDIs. The modified model correctly simulates (i) a higher value of Cmax on day 1 for higher biovalability; (ii) a higher PRED value at steady state for lower apparent clearance effect CLSCALE/FSCALE values.

Dose delay

Deviations in the administration schedule of a drug are simulated including three parameters in the dataset: the dose number DOSCUMN, obtained with NMdata::addTAPD(), the specific number of the delayed dose DELAYDOS and the time delay ALAG.

dat.sim.alag = data.sim |> NMexpandDoses() |> addTAPD() # expand doses and add dose number
dat.sim.alag[,ROW:=.I] # restore ROW with correct value
# add delayed dose
dat.sim.dos.delay=dat.sim.alag[,data.table(DELAYDOS=c(2,3,4,5,6,7))
                                ,by=dat.sim.alag]
# add dose delay 
dat.sim.alag.final=dat.sim.dos.delay[,data.table(ALAG=c(0,6,12,18,24))
                                     ,by=dat.sim.dos.delay]
# update ID 
dat.sim.alag.final[,ID:=.GRP,by=.(ALAG,DELAYDOS,ID)]
setorder(dat.sim.alag.final,ID,TIME,EVID)
#dat.sim.alag.final = dat.sim.alag.final |> fill(DOSCUMN)

The following patient has a 6 hours delay to the administration of dose 2.

TIME AMT DOSCUMN DELAYDOS ALAG
0 300 1 2 6
24 150 2 2 6
48 150 3 2 6
72 150 4 2 6
96 150 5 2 6
120 150 6 2 6
144 150 7 2 6

The time delay is included in the modified control stream with NMsim adding a single line at the end of the PK section, where the dose delay on compartment 1 ALAG1 is modified for the dose with dose number DOSCUMN equal to the target delayed dose number DELAYDOS

simres.alag <- NMsim(file.mod=file.mod
              ,data=dat.sim.alag.final
              ,name.sim="alag"
              ,modify.model=list(PK=
                    add("IF(DOSCUMN.EQ.DELAYDOS) ALAG1=ALAG")))

The effect on concentration-time profiles and daily AUC are shown in the two images below.

Effect of time delay (0, 12 and 24 hours on dose 2) on concentration (`PRED`) profile as a function of time computed by `NMsim` modified model. The implementation of the modified model simply consists of the addition of the variables `DOSCUMN`, `DELAYDOS`, and `ALAG` to the original data set, and the addition of one line of code to the PK section of the control stream via `modify.model`.

Effect of time delay (0, 12 and 24 hours on dose 2) on concentration (PRED) profile as a function of time computed by NMsim modified model. The implementation of the modified model simply consists of the addition of the variables DOSCUMN, DELAYDOS, and ALAG to the original data set, and the addition of one line of code to the PK section of the control stream via modify.model.

Daily exposure on day 3 as a function of time delay for dose 2 (left) and dose 3 (right). The simulation results predict an increased risk for possible safety concerns (left panel, over-exposure) and loss of efficacy (right panel, under-exposure) as dose time delay gets larger.

Daily exposure on day 3 as a function of time delay for dose 2 (left) and dose 3 (right). The simulation results predict an increased risk for possible safety concerns (left panel, over-exposure) and loss of efficacy (right panel, under-exposure) as dose time delay gets larger.

AUC

The following example computes daily exposure:

  • at post-processing, using trapezoidal method, after the unmodified NONMEM model is simulated on three different fine grids with evenly spaced time steps of 0.25, 1, and 4 hours, respectively, labelled AUC trapez 
file.mod.auc <- system.file("examples/nonmem/xgxr046.mod",
                        package="NMsim")
data.sim.auc <- NMreadCsv(system.file("examples/derived/dat_sim1.csv",
                                 package="NMsim"))
data.sim.auc[,AMT:=1000*AMT]

# time step 1hr
data.sim.1hr=data.sim.auc
data.sim.1hr[,TSTEP:="1hr"]

# time step 0.25hr
data.sim.0.25hr=addEVID2(data.sim.auc[EVID==1],CMT=2,time.sim=seq(0,192,by=0.25))
data.sim.0.25hr[,TSTEP:="0.25hr"]

# time step 4hr
data.sim.4hr <- addEVID2(data.sim.auc[EVID==1],CMT=2,time.sim=seq(0,192,by=4))
data.sim.4hr[,TSTEP:="4hr"]

sres.trapez <- NMsim(file.mod=file.mod.auc
      ,data=list(data.sim.1hr # run NMsim on a list of data sets to run all different scenarios at once
                 ,data.sim.0.25hr
                 ,data.sim.4hr)
      ,seed=12345
      ,table.vars=cc(PRED,IPRED)
      ,name.sim="AUC.trapez"
      ,reuse.results=reuse.results
      )

# daily AUC computation
sim.auc=sres.trapez[EVID==2,.(ID,TIME,PRED,time.step=TSTEP)] 
sim.auc[,DAY:=(TIME%/%24)+1] # define DAY variable
# create duplicate time steps for end-of-day
# (e.g 24h belongs to both day 1 and day2)
sim.dupli.24h=sim.auc[TIME%%24==0 & TIME>0,.(ID
                                             ,DAY=DAY-1
                                             ,TIME
                                             ,PRED
                                             ,time.step)]
sim.auc.final <- rbind(sim.auc
                       ,sim.dupli.24h) |> setorder(ID,DAY,TIME,time.step)

sim.auc.trapez<-sim.auc.final[DAY<8,.(AUC=NMcalc::trapez(TIME,PRED))
                             ,by=.(ID,DAY,time.step)]
# stmp=sres1[EVID==2,.(ID,TIME,PRED)]
# stmp[,DAY:=(TIME%/%24)+1]
# sAUC<-stmp[,.(AUC=NMcalc::trapez(TIME,PRED)),by=.(ID,DAY)]
  • at run time, with an NMsim modified script, on a course grid with evenly spaced time steps of 24 hours, labelled AUC $DES. Of note, this task requires additions to the control stream in multiple sections.
# AUC with NMsim -  time step 24hr
data.sim2 <- addEVID2(data.sim.auc[EVID==1],CMT=2,time.sim=seq(0,by=24,length.out=9))

sres.des <- NMsim(file.mod=file.mod.auc
          ,data=data.sim2
          ,table.vars=cc(PRED,IPRED,AUCNMSIM)
          ,name.sim="AUC.nmsim"
          ,seed=12345
          ,reuse.results=reuse.results
          ,modify.model=list(MODEL=add("COMP=(AUC)")
                  ,DES=add("DADT(3)=A(2)/V2")
                  ,ERROR=add("AUCCUM=A(3)"
                             ,"IF(NEWIND.NE.2) OLDAUCCUM=0"
                             ,"AUCNMSIM = AUCCUM-OLDAUCCUM"
                             ,"OLDAUCCUM = AUCCUM")))

sres.des[,DAY:=TIME%/%24]
sres.des[,AUC.NMsim:=AUCNMSIM]
sres.auc.des=sres.des[AUCNMSIM!=0,.(TIME,DAY,PRED,AUC.NMsim)]

The plot comparing the DES and trapez AUC is obtained with the code below

sres.final=mergeCheck(sim.auc.trapez[,.(DAY,AUC,time.step)]
                      ,sres.auc.des[,.(DAY,AUC.NMsim)]
                      ,by="DAY")

ggplot(data=sres.final,aes(AUC.NMsim,AUC,colour=factor(time.step)))+
  geom_point(size=4)+
  labs(colour="Time step in fine grid")+
  scale_color_manual(values=c("orange", "blue", "darkgreen"))+
  xlim(c(0,17))+
  ylim(c(0,17))+
  ylab("AUC 0-24h by trapez method, on fine grid")+
  xlab("AUC 0-24h by integration through $DES, on coarse grid")+
  geom_abline(slope=1, intercept=0)
Daily exposures computed at run time ($DES, coarse grid, x-axis) and post-processing time (`trapez`, fine grids, y-axis). AUC (trapez) converges to the value computed with $DES method as the time step is reduced. Includes identity line.

Daily exposures computed at run time ($DES, coarse grid, x-axis) and post-processing time (trapez, fine grids, y-axis). AUC (trapez) converges to the value computed with $DES method as the time step is reduced. Includes identity line.