Data Source

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: Simulated_gentamicin_pk.csv }

Data Preprocessing

Because PoPy expects data in a slightly different format to Nonmem, we make some modifications to the incoming data using the PREPROCESS block of the input script.

The critical columns for the original Nonmem data set looks something like Table 55.

and are missing a TYPE column along with a flag column that serves essentially the same purpose as Nonmem’s MDV column.

We add these missing columns using the PREPROCESS section of the fit to modify every row of the incoming data file as it is read:

PREPROCESS: |
    # Call this helper function that (a) converts Nonmem's EVID column
    # values into PoPy's TYPE column values and (b) adds a DV_FLAG
    # column that is 1 for observation rows (EVID=0) and 0 everywhere else.
    nonmem2popy()

    # Augment PoPy's new TYPE column value for dose rows (where AMT>0)
    # with information about whether the dose is an infusion (RATE>0)
    # or a bolus (RATE=0).
    if c[AMT] > 0:
        if c[RATE] > 0:
            c[TYPE] += ":DOSE_infrate"
        else:
            c[TYPE] += ":DOSE_bolus"

    # Make a copy of Nonmem's DV column to give it a more descriptive
    # name, DRUG_CONC, and update the corresponding flag.
    # (This is optional and could be left out after changing the
    # PREDICTIONS block to use DV rather than DRUG_CONC.)
    c[DRUG_CONC] = c[DV]
    c[DRUG_CONC_FLAG] = c[DV_FLAG]

Technically the Nonmem data set above has no CMT field, but CMT defaults to 1 in Nonmem if it is not provided. PoPy does not assign a compartment in the data file, preferring to assign the compartment in the model definition (where compartments have meaning).

The (critical columns of the) preprocessed dataset used by PoPy will now look like that in Table 56.

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 Nonmem’s’ MDV flag. The new TYPE field is formatted to indicate either bolus or infusion doses in PoPy.

Fixed/Random Effects Specification

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 been 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 ``SIGMA``s are treated as fixed effects and also defined in EFFECTS.

PK Model Parameters Specification

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 PoPy’s’ 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.

ODEs Specification

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.

Likelihood Error Specification

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 (even though it is not clear - to the PoPy developers, at least - 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.

Parameter Fitting Methods Specification

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.