Overview¶

A real world model, based on a gentamicin PK study [Harling2015]. It has a combined bolus + infusion dose and the original PK model uses ADVAN3 TRANS4 in Nonmem. There are some covariate effects on the clearance and volume of distribution of the central compartment parameters. There are 210 individuals and 1949 data rows. See DDMoRe: 0061

Data Conversion¶

We describe the conversion of this data set in Nonmem format:-

c:\PoPy\validation\ddmore0061\Simulated_gentamicin_pk.csv

To this data set in PoPy format:-

c:\PoPy\validation\ddmore0061\popy\popy_data.csv

Using the N2PDat Script located here:-

c:\PoPy\validation\ddmore0061\popy\n2p_script.pyml

The original Nonmem data set looks something like Table 60.

The critical entries for the N2PDat Script are as follows:-

INPUT_NONMEM_FIELDS:
     obs_cmt_numbers: [1]
     dose_cmt_numbers: [1]

These two entries specify that the observations and doses all occur in compartment one. Technically the Nonmem data set above has no CMT field. However CMT defaults to one in Nonmem if it is NOT provided.

The dataset output by N2PDat Script is the same as the input, but with the 3 extra columns shown in Table 61.

Here the new ‘DRUG_CONC’ field will be the target value for the new PoPy script. The ‘DRUG_CONC_FLAG’ switches the observations on/off appropriately. Note that ‘DRUG_CONC_FLAG’ is the inverse of the Nonmem ‘MDV’ flag. The new ‘TYPE’ field is formatted to process either bolus or infusion doses in PoPy.

Script Conversion¶

We describe the manual conversion of the Nonmem script file:-

c:\PoPy\validation\ddmore0061\Executable_gentamicin_pk.mod

To create the PoPy Fit Script here:-

c:\PoPy\validation\ddmore0061\popy\fit_script.pyml

The Nonmem script loads the data as follows:-

$DATA Simulated_gentamicin_pk.csv IGNORE=@

Whereas the PoPy script uses this syntax:-

FILE_PATHS: {input_data_file: popy_data.csv}

In the Nonmem script, the theta/omega/sigma fixed effect definitions are:-

$THETA
    (0,4.1) ; CL
    (1,8.63) ; V1
    (0,1.21) ; Q
    (0,8.12) ; V2

    (-0.01,0.00982,0.016) ; BCLC
    (-0.409) ; ALB
    (0,0.00683) ; DCLC

$OMEGA
    0.0465 ; CL
    0.376 ; Q

$SIGMA
    0.0297 ; E1
    0.0486 ; E2

In PoPy these fixed effect are defined in EFFECTS:-

EFFECTS:
    POP: |
        f[CL] ~ P 4.1
        f[V1] ~ unif(1, +inf) 8.63
        f[Q] ~ P 1.21
        f[V2] ~ P 8.12
        f[BCLC_eff] ~ unif(-0.01, 0.016) 0.00982
        f[ALB_eff] = -0.409
        f[DCLC_eff] ~ P 0.00683

        f[CL_isv, Q_isv] ~ diag_matrix() [ [ 0.0465, 0.376 ] ]

        f[PNOISE_var] ~ P 0.0297
        f[ANOISE_var] ~ P 0.0486

Note the ranges and initial values of the f[X] variables are defined using different syntax from the $THETA variables. Note also that the f[X] have explicit names, whereas the Nonmem list of thetas has to be clumsily annotated with comments to make it human readable.

In PoPy you can specify a diagonal matrix, which can be explicitly used as a covariance matrix input to a Multivariate Normal Distribution distributed r[X] vector, as follows:-

EFFECTS:
    ID: |
        r[CL_isv, Q_isv] ~ mnorm([0,0], f[CL_isv, Q_isv])

There r[CL_isv] is the equivalent of ‘ETA(1)’ and r[Q_isv] is the equivalent of ‘ETA(2)’. In Nonmem you just have to remember the index values, whereas PoPy uses actual real life variable names.

In PoPy the sigmas are fixed effects and also defined in EFFECTS.

The Nonmem PK section is as follows:-

$PK
    BSA = 71.84*(WT**0.425)*(HT**0.725)
    BSA = BSA/10000
    HT1 = HT/100
    BMI = WT/HT1**2

    IF(NEWIND.NE.2) BCLC = CLC
    DCLC = CLC-BCLC
    IF(NEWIND.NE.2) BBSA = BSA

    TVCL = THETA(1) * (1 + THETA(5)*(BCLC-81) + THETA(7)*DCLC)
    TVV1 = THETA(2) * BBSA*(ALB/34)**THETA(6)
    TVQ  = THETA(3)
    TVV2 = THETA(4)

    CL = TVCL * EXP(ETA(1))
    V1 = TVV1
    Q  = TVQ  * EXP(ETA(2))
    V2 = TVV2

    S1 = V1

This maps almost exactly to the PoPy MODEL_PARAMS section:-

MODEL_PARAMS: |
    BSA = 71.84 * (c[WT]**0.425) * (c[HT]**0.725)
    BSA = BSA / 10000
    HT1 = c[HT] / 100
    BMI = c[WT] / HT1**2

    BCLC = c[CLC]
    DCLC = c[CLC] - BCLC
    BBSA = BSA

    TVCL = f[CL] * (1 + f[BCLC_eff]*(BCLC-81) + f[DCLC_eff]*DCLC)
    TVV1 = f[V1] * BBSA * (c[ALB]/34)**f[ALB_eff]
    TVQ  = f[Q]
    TVV2 = f[V2]

    m[CL] = TVCL * EXP(r[CL_isv])
    m[V1] = TVV1
    m[Q]  = TVQ  * EXP(r[Q_isv])
    m[V2] = TVV2

    m[PNOISE_var] = f[PNOISE_var]
    m[ANOISE_var] = f[ANOISE_var]

Note in the above Nonmem specifies ‘S1 = V1’, to scale the amounts in compartment one. This is not necessary in PoPy, as the compartment amounts are converted to concentrations in the PREDICTIONS section explicitly.

Also note that PoPy requires the lines:-

m[PNOISE_var] = f[PNOISE_var]
m[ANOISE_var] = f[ANOISE_var]

To propagate the f[X] noise fixed effects which are similar to the Nonmem builtin sigma variables.

Note in this example the Nonmem ‘NEWIND’ keyword is superfluous, as the c[CLC] covariate and BSA variable are constant within each individual.

The Nonmem control file prescribes the compartment model on this line:-

$SUBROUTINE ADVAN3 TRANS4

That uses the inbuilt ‘ADVAN3’ model with ‘TRANS4’ parametrisations. This uses the analytic solution for a two compartment model that expects the CL, V1, Q, V2 parameters to be defined in the Nonmem PK section.

The PoPy control file here specifies the compartment model equations explicitly in the DERIVATIVES section:-

DERIVATIVES: |
    dose[DOSE_bolus] = @bolus{amt:c[AMT], lag:0.0}
    dose[DOSE_infrate] = @inf_rate{amt:c[AMT], lag:0.0, rate:c[RATE]}
    d[CENTRAL] = dose[DOSE_bolus] + dose[DOSE_infrate] - s[CENTRAL]*m[Q]/m[V1] + s[PERI]*m[Q]/m[V2] - s[CENTRAL]*m[CL]/m[V1]
    d[PERI]    =                                         s[CENTRAL]*m[Q]/m[V1] - s[PERI]*m[Q]/m[V2]

Note that PoPy has an equivalent builtin version of ‘ADVAN3 TRANS4’ as follows:-

DERIVATIVES: |
    s[CENTRAL, PERI] = @iv_two_cmp_cl{ dose: @bolus{amt:c[AMT]}, CL: m[CL], V1: m[V1], Q: m[Q], V2: m[V2] }

However this function is unable to process two different types of dose (bolus + infusion), unlike the long hand version above. Note that PoPy will be using an ordinary differential equation solver here (instead of a built in analytic solution), so we also need to specify this PoPy field:-

ODE_SOLVER: {CPPLSODA: {}}

Which tells PoPy to use the c++ LSODA solver, the same method as Nonmem Advan13.

The Nonmem ERROR section is as follows:-

$ERROR
    Y = F*EXP(ERR(1)) + ERR(2)
    IPRED=F
    IRES=DV-IPRED

The PoPy equivalent is shown below in the PREDICTIONS section:-

PREDICTIONS: |
    p[DV_CENTRAL] = s[CENTRAL]/m[V1]
    var = m[PNOISE_var] * p[DV_CENTRAL]**2 + m[ANOISE_var]
    c[DRUG_CONC] ~ norm(p[DV_CENTRAL], var)

Note here the line ‘Y = F*EXP(ERR(1)) + ERR(2)’ defines an proportional exponential normal noise model with an additional additive normal noise component. It is not clear (to the PoPy developers) what the distribution of a log normal + a normal distribution is. PoPy requires a known distribution to calculate a likelihood. If you approximate the exponential (assuming ERR(1) is small) and get ‘Y = F*(1+ ERR(1)) + ERR(2)’. Then this is the standard proportional + additive noise model. In this case, the PoPy Fit Script return a very similar objective function to the DDMoRe Nonmem objective function. Indicating that Nonmem also makes the same approximation.

The Nonmem control file specifies using FOCE fitting (METHOD=1) below:-

$ESTIMATION MAXEVALS=0 SIG=3 PRINT=10 METHOD=1 INTER

The PoPy equivalent, specifying the ND fitting method is:-

FIT_METHODS: [ND: {max_n_main_iterations: 0}]

Setting ‘max_n_main_iterations’ to zero, means that PoPy will optimise the r[X] parameters and leave the f[X] unchanged.