NONMEM Scaling factor 分布体积编号问题 Weibull吸收模型示例

在使用 NONMEM 时注意确保scaling factor对应体积的正确 （尤其当ADVAN 6等自定义模型时）

记录一次NONMEM bug shooting：

场景： 复写文献里的一个 ** 2房室，Weibull吸收，一级代谢的 imatinib popPK。**（模型来自Park, G. J., Park, W. S., Bae, S., Park, S. M., Han, S., & Yim, D. S. (2016). Population pharmacokinetics of imatinib mesylate in healthy Korean subjects. Translational and Clinical Pharmacology, 24(2), 96–104.）
模型结构：

模型参数 (去掉Age covariate）：

parame	值
CL (L/h)	13.6
Vc (L)	153
Vp (L)	64
Q (L/h)	8.64
K1	0.998
GAMMA	2.24

模拟数据集中DOSE是mg，在无scaling factor情况下模型得到的浓度为 mg/L, 1mg/L=1000ng/mL

代码1 （错误示范）

细心的小伙伴一看代码就知道哪里错了

$PROB ;;2-comp Weibull abosorption
$DATA imatinib_sim_children_SS2_form_vpcPLOT.csv IGNORE=@
$INPUT ID,TIME,AMT,RATE=DROP,DV,ROUTE=DROP,OCC=DROP,MDV,EVID,CMT=DROP,WT,AGE=DROP,SEX,HT,II,SS,TAD
$SUBROUTINE ADVAN6 TRANS1 TOL=6
$MODEL 
COMP(DEPOT,DEFDOS)
COMP(CENT, DEFOBS)
COMP(PERIP)

$PK
TVCL     = THETA(1)
TVV1     = THETA(2)
TVV2     = THETA(3)
TVQ      = THETA(4)
TVKA1    = THETA(5)
TVGAMMAP = THETA(6)

CL     = TVCL*EXP(ETA(1))
V1     = TVV1*EXP(ETA(2))
V2     = TVV2*EXP(ETA(3))
Q      = TVQ*EXP(ETA(4))
KA1    = TVKA1*EXP(ETA(5))
GAMMAP = TVGAMMAP*EXP(ETA(6))

K12  = 1 - EXP(-(KA1*TAD)**GAMMAP)
K20  = CL/V1
K23  = Q/V1
K32  = Q/V2

S2 = V2/1000

$DES
DADT(1) = -K12*A(1)
DADT(2) = K12*A(1)+K32*A(3)-K20*A(2)-K23*A(2)
DADT(3) = K23*A(2)-K32*A(3)

$ERROR
IPRED = F
Y     = IPRED*(1+EPS(1))+EPS(2)
IRES  = DV - IPRED
SD    = SQRT(SIGMA(1,1)*(IPRED)**2 + SIGMA(2,2))
IF(SD.EQ.0)  SD=1
IWRES = IRES/SD

$THETA 13.6  FIX ;CL     1
$THETA 153   FIX ;V1     2
$THETA 64    FIX ;V2     3
$THETA 8.64  FIX ;Q      4
$THETA 0.998 FIX ;KA1    5
$THETA 2.24  FIX ;GAMMAP 6

$OMEGA 0 FIX
$OMEGA 0 FIX
$OMEGA 0 FIX
$OMEGA 0 FIX
$OMEGA 0 FIX
$OMEGA 0 FIX

$SIGMA 0.108 ;PROP ERROR FOR PK
$SIGMA 0.0001;ADD ERROR FOR PK

$SIMULATION (666)  ONLYSIM SUBPROBLEMS=300
$TABLE ID TIME TAD II SS AMT DV MDV EVID WT AGE SEX HT CL TVV1 V1 V2 Q KA1 GAMMAP IPRED Y IRES SD IWRES EP1 EP2 NOPRINT ONEHEADER FILE=Park_2016_simulate_population_percentiles.fit

此时得到的percentile plot 浓度高到非常离谱：
（从上向下的线依次为模拟300个datasets在各时间点浓度数据的 0.95, 0.75, 0.5, 0.25, 0.05分位数。）

代码2 （正确示范）

$PROB ;;2-comp Weibull abosorption
$DATA imatinib_sim_children_SS2_form_vpcPLOT.csv IGNORE=@
$INPUT ID,TIME,AMT,RATE=DROP,DV,ROUTE=DROP,OCC=DROP,MDV,EVID,CMT=DROP,WT,AGE=DROP,SEX,HT,II,SS,TAD
$SUBROUTINE ADVAN6 TRANS1 TOL=6
$MODEL 
COMP(DEPOT,DEFDOS)
COMP(CENT, DEFOBS)
COMP(PERIP)

$PK
TVCL     = THETA(1)
TVV1     = THETA(2)
TVV2     = THETA(3)
TVQ      = THETA(4)
TVKA1    = THETA(5)
TVGAMMAP = THETA(6)

CL     = TVCL*EXP(ETA(1))
V1     = TVV1*EXP(ETA(2))
V2     = TVV2*EXP(ETA(3))
Q      = TVQ*EXP(ETA(4))
KA1    = TVKA1*EXP(ETA(5))
GAMMAP = TVGAMMAP*EXP(ETA(6))

K12  = 1 - EXP(-(KA1*TAD)**GAMMAP)
K20  = CL/V1
K23  = Q/V1
K32  = Q/V2


$DES
DADT(1) = -K12*A(1)
DADT(2) = K12*A(1)+K32*A(3)-K20*A(2)-K23*A(2)
DADT(3) = K23*A(2)-K32*A(3)

$ERROR
IPRED = A(2)/V1
Y     = IPRED*(1+EPS(1))+EPS(2)
IRES  = DV - IPRED
SD    = SQRT(SIGMA(1,1)*(IPRED)**2 + SIGMA(2,2))
IF(SD.EQ.0)  SD=1
IWRES = IRES/SD

$THETA 13.6  FIX ;CL     1
$THETA 153   FIX ;V1     2
$THETA 64    FIX ;V2     3
$THETA 8.64  FIX ;Q      4
$THETA 0.998 FIX ;KA1    5
$THETA 2.24  FIX ;GAMMAP 6

$OMEGA 0 FIX
$OMEGA 0 FIX
$OMEGA 0 FIX
$OMEGA 0 FIX
$OMEGA 0 FIX
$OMEGA 0 FIX

$SIGMA 0.108 ;PROP ERROR FOR PK
$SIGMA 0.0001;ADD ERROR FOR PK

$SIMULATION (666)  ONLYSIM SUBPROBLEMS=300

此时我们没有用 scaling factor S2，得到的浓度为mg/L，乘1000后得到ng/mL。 percentile plot:

原因分析

在代码1的模型error 中插入以下代码，并输出table

$ERROR
EP1 = EPS(1)
EP2 = EPS(2)
IPF = F
IPA1 = A(1)
IPA2 = A(2)
IPA1V = A(1)/V1
IPA2V = A(2)/V2
$TABLE ID IPF IPA1 IPA2 IPA1V IPA2V IPRED Y MDV TAD EP1 EP2 NOPRINT ONEHEADER FILE=check.fit

See check.fit (head)

可以看到事实上我们的F(IPF) 等于了 IPAV*1000，即在代码1中 F = A(2)/V2。A(2)是我们定义的defobs 观测室，但V2是我们给的周边室的分布体积Vp (我们$MODEL 中定义的第三个COMP)。 我们想要的正确的观测 应该是F = A(2)/V1，即代码2中的IPRED。

总结

我们知道当我们不写明scaling factor时，F会输出观测室的物质的量（即代码1中如果没有S2=V2/1000，F的输出会是A(2)）。在nonmem中我们定义了三个房室，分布体积的参数只有V1和V2，即事实上V1是我们定义的nonmem中第二个房室（中央室）的分布体积，V2是我们定义的的nonmem中第3个房室（周边室）的分布体积。在写代码的时候经常容易搞错，以至于长时间检查无法发现，个人建议在写nonmem时[1]可以尽量用VC，VP（VP1,VP2) 而非V1，V2, V3。[2] 尤其在使用Advan6 自定义模型时可以尽量使用IPRED 直接定义而少用F。

Weibull 吸收模型

上述的代码2可以是一个Weibull吸收模型的写法示范，percentile图相对浓度较高是因为图中描述的是在稳态给药后的血药浓度。下图是一个第一天给药的血药浓度percentile plot

对应参考Berkeley Madonna的结果

# Berkeley Madonna 代码
METHOD RK4
STARTTIME = 0
STOPTIME = 100
DT = 0.02

{Type Equations Here.}
Dose = 400
V1 = 153
CL = 13.6
index = (0.998*time)^2.24
ka =1-EXP(-index)
V2 = 64
Q = 8.64

k10 = CL/V1
k12 = Q/V1
k21 = Q/V2

d/dt(D) = -ka*D
d/dt(A1) = ka*D-k10*A1- k12*A1+k21*A2 
d/dt(A2) = k12*A1 - k21*A2

init D = Dose
init A1 = 0
init A2 = 0

CP = A1/V1