Overview¶

A real world model, based on a QT study [Cheung2015]. The model incorporates a circadian function with no compartment derivatives. There are 59 individuals and 5790 total data rows. See DDMoRe: 0093

Data Conversion¶

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

c:\PoPy\validation\ddmore0093\Simulated_dataset.csv

To this data set in PoPy format:-

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

Using the N2PDat Script located here:-

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

The final data set looks something like Table 62.

Note the main output of the data conversion is the ‘TYPE’ field which has the correct entries for PoPy corresponding to the Nonmem ‘EVID’ field. as follows:-

2 -> pred
0 -> obs
1 -> dose:dose_bolus

Note however that in this model the dosing is not actually used in the Nonmem control file.

Script Conversion¶

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

c:\PoPy\validation\ddmore0093\Executable_run7b.ctl

To create the PoPy Fit Script here:-

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

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.

Whereas the PoPy script uses this syntax:-

FILE_PATHS: {input_data_file: output_popy_data.csv}
PREPROCESS: |
    if c[TYPE] == 'dose:dose_bolus': return

Which loads the PoPy data and also excludes the dosing rows.

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: |
        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
        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 ]
        ]
        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:-

$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.

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.

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}]