Data Source & Preprocessing

The Nonmem script loads the data as follows:-

$DATA Simulated_dataset.csv
   IGNORE = #
   IGNORE = (EVID == 1)

Note the ‘IGNORE’ syntax removes the dosing rows from the data set, effectively preprocessing the dataset at the point of loading.

We achieve the same thing in PoPy using the PREPROCESS block:

FILE_PATHS: { input_data_file: Simulated_dataset.csv }
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()

    # Suppress warning about having doses we don't use
    # by ignoring dose rows (EVID==1) altogether.
    if c[TYPE] == 'dose': return

which loads the PoPy data, converts EVID values to TYPE values and also excludes the dosing rows.

Fixed/Random Effects Specification

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

$THETA (0, 372.6)      ; 1 Baseline QTcF [ms]
$THETA (0, 0.01844)    ; 2 Amplitude  24h
$THETA (0, 3.62)       ; 3 Peak shift 24h
$THETA (0, 0.01392)    ; 4 Amplitude  12h circadian rhythm
$THETA (0, 1.301)      ; 5 Peak shift 12h circadian rhythm
$THETA (0, 0.003441)   ; 6 Slope (linear effect) parameter

$OMEGA  0.0392         ; 1 Baseline QTcF - (BSV)
$OMEGA  0.3523         ; 2 Amplitude 24h  - BSV
$OMEGA  2.007          ; 3 Peak shift 24h - BSV
$OMEGA  0.3848         ; 4 Amplitude 12h  - BSV
$OMEGA  0.993          ; 5 Peak shift 12h - BSV
$OMEGA  0.0001         ; 6 Slope          - BSV

$SIGMA  0.01802        ; 1 Proportional residual error [ms]

The equivalent f[X] in the PoPy script are:-

EFFECTS:
    POP: |
        # THETAs
        f[BASE] ~ P 372.6
        f[AMP24] ~ P 0.01844
        f[SHFT24] ~ P 3.62
        f[AMP12] ~ P 0.01392
        f[SHFT12] ~ P 1.301
        f[SLOPE] ~ P 0.003441

        # OMEGAs
        f[BASE_bsv,AMP24_bsv,SHFT24_bsv,AMP12_bsv,SHFT12_bsv,SLOPE_bsv] ~ diag_matrix() [
            [ 0.0392, 0.3523, 2.007, 0.3848, 0.993, 0.0001 ]
        ]

        # SIGMAs
        f[NOISEVAR] ~ P 0.01802
    ...

Additionally PoPy has the ‘ID’ section of the EFFECTS:-

EFFECTS:
    ...
    ID: |
        r[ BASE, AMP24, SHFT24, AMP12, SHFT12, SLOPE ] ~ mnorm(
            [0,0,0,0,0,0],
            f[BASE_bsv,AMP24_bsv,SHFT24_bsv,AMP12_bsv,SHFT12_bsv,SLOPE_bsv]
        )

This section defines r[BASE] and r[AMP24] etc. These r[X] definitions are a more explicit version of Nonmem implicitly creating ETA(X) variables for each of the omega variables. With Nonmem you have to remember what each of the ETA indices represent. Hence the heavily commented ‘PRED’ block below, which defined the variables for each individual in the Nonmem script.

PK Model Parameters Specification

$PRED
   BASE   = THETA(1) * EXP(ETA(1))                       ; Baseline QTcF with between subject varaibility (BSV)
   AMP24  = THETA(2) * EXP(ETA(2))                       ; The amplitude of the 24h circadian rhythm
   SHFT24 = THETA(3) + ETA(3)                            ; The peak shift of the 24h circadian rhythm
   AMP12  = THETA(4) * EXP(ETA(4))                       ; The amplitude of the 12h circadian rhythm
   SHFT12 = THETA(5) + ETA(5)                            ; The peak shift of the 12h circadian rhythm

   CIRC24 = AMP24 * COS(2 * 3.14 * (TIME - SHFT24)/24)   ; 24 hour circadian rhythm
   CIRC12 = AMP12 * COS(2 * 3.14 * (TIME - SHFT12)/12)   ; 12 hour circadian rhythm

   RYTM   = BASE * (1 + CIRC24 + CIRC12)                 ; Change in baseline QTcF over the day due to circadian rhythm

   SLOPE  = THETA(6) + ETA(6)                            ; Linear effect
   EFF    = SLOPE * CPP                                  ; Linear effect

   ...

The MODEL_PARAMS section of the PoPy script is very similar:-

MODEL_PARAMS: |

    BASE   = f[BASE] * exp(r[BASE])
    AMP24  = f[AMP24] * exp(r[AMP24])
    SHFT24 = f[SHFT24] + r[SHFT24]
    AMP12  = f[AMP12] * exp(r[AMP12])
    SHFT12 = f[SHFT12] + r[SHFT12]

    CIRC24 = AMP24 * COS(2 * 3.14 * (c[TIME] - SHFT24)/24)
    CIRC12 = AMP12 * COS(2 * 3.14 * (c[TIME] - SHFT12)/12)
    m[RYTM] = BASE * (1 + CIRC24 + CIRC12)

    SLOPE = f[SLOPE] + r[SLOPE]
    m[EFF] = SLOPE * c[CPP]
    m[NOISEVAR] = f[NOISEVAR]

The above is almost identical apart from the line ‘m[NOISEVAR] = f[NOISEVAR]’. Note in MODEL_PARAMS it is necessary to define m[X] variables that you wish to use in subsequent sections, e.g. DERIVATIVES or PREDICTIONS.

ODEs Specification

This model does not require or use differential equations so the $DES section is absent in the Nonmem file and the DERIVATIVES section is absent from the PoPy file.

Likelihood Error Specification

The Nonmem error model (also here part of the ‘PRED’ section) is defined as follows:-

IPRED  = RYTM  + EFF        ; Linear direct effect model
W      = IPRED * SIGMA(1,1)
IWRES  = (QTCF - IPRED)/W

Y      = IPRED + IPRED*EPS(1)

The equivalent in PoPy is the PREDICTIONS block:-

PREDICTIONS: |
    p[CONC_PLASMA] = m[RYTM] + m[EFF]
    var = m[NOISEVAR] * p[CONC_PLASMA]**2
    c[QTCF] ~ norm(p[CONC_PLASMA], var)

Note PoPy explicitly defines the likelihood by comparing c[QTCF] from the data file with the predicition p[CONC_PLASMA], using a normal distribution with proportional variance.

This is more explicit than the Nonmem line ‘Y = IPRED + IPRED*EPS(1)’, which implicitly uses the ‘DV’ Nonmem magic variable for Y. The likelihood is expressed as a function of EPS normal distributions.

Parameter Fitting Methods Specification

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

$ESTIMATION METHOD=1 MAXEVALS=99999 INTER NOABORT PRINT=5

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

FIT_METHODS: [ND: {max_n_main_iterations: 100}]