Overview¶

A real world model, based on a gentamicin PK study [Germovsek2017]. A Population PK model with IOV, infusion dosing and custom derivative equations to implement the PK. The dataset consists of 205 individuals and 2788 rows. See DDMoRe: 0238

Data Conversion¶

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

c:\PoPy\validation\ddmore0238\Simulated_simdataDDM.csv

To this data set in PoPy format:-

c:\PoPy\validation\ddmore0238\popy\output_popy_data.csv

Using the N2PDat Script located here:-

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

The original Nonmem data set looks something like Table 63.

The final data set looks something like Table 64.

Here the PoPy ‘TYPE’ field is a more explicit version of the Nonmem ‘EVID’ field. The PoPy ‘DV_CENTRAL’ field is a copy of the Nonmem ‘DV’ field, with an additional ‘DV_CENTRAL_FLAG’ field to make the true observed values more obvious. Note here PoPy could also use the Nonmem ‘DV’ field because it is always defined for all rows with c[EVID] =0 or equivalently c[TYPE] =’obs’.

Script Conversion¶

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

c:\PoPy\validation\ddmore0238\Executable_run35b_ddm2.ctl

To create the PoPy Fit Script here:-

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

The Nonmem script loads the data as follows:-

$DATA simdataDDM.csv IGNORE=@

Whereas the PoPy script uses this syntax:-

FILE_PATHS:
    # path to input comma separated value file in popy data format
    input_data_file: output_popy_data.csv

In the Nonmem script, the theta variable definitions are:-

$THETA  (0,6.20684) ; 1. TVCL (lower bound,initial estimate)
$THETA  (0,26.5004) ; 2. TVV1  (lower bound,initial estimate)
$THETA  (0,2.15099) ; 3. TVQ
$THETA  (0,21.151) ; 4. TVV2
$THETA  (0,0.270697) ; 5. TVQ2
$THETA  (0,147.893) ; 6. TVV3
$THETA  55.4 FIX ; 7. T50
$THETA  3.33 FIX ; 8. Hill
$THETA  -0.129934 ; 9. power exponent on creatinine
$THETA  (0,1.70302) ; 10. PNA50

In the PoPy script the equivalent f[X] variables are:-

EFFECTS:
    POP: |
        f[TVCL] ~ P 6.20684
        f[TVV1] ~ P 26.5004
        f[TVQ] ~ P 2.15099
        f[TVV2] ~ P 21.151
        f[TVQ2] ~ P 0.270697
        f[TVV3] ~ P 147.893
        f[T50] = 55.4
        f[HILL] = 3.33
        f[CRPWR] = -0.129934
        f[P50] ~ P 1.70302

Notice that PoPy uses the natural equal sign for f[X] that are fixed and the ‘~’ for f[X] that need to be estimated.

In the Nonmem script, the omega variable definitions are:-

$OMEGA  BLOCK(2)
 0.175278  ; variance for ETA(1), initial estimate
 0.115896 0.112362  ; COvariance ETA(1)-ETA(2), var for ETA(2), initial estimate
$OMEGA  0  FIX
$OMEGA  0.131759
$OMEGA  0  FIX
$OMEGA  0.177214
$OMEGA  BLOCK(1)
 0.0140684  ;  7. IOV_CL
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME

The ugly code above defines a 2x2 matrix of omega fixed effects, followed by 4 scalar omega fixed effects. The ‘$OMEGA BLOCK(1)’ + ‘$OMEGA BLOCK(1) SAME’ lines defines a single IOV omega fixed effects to be estimated over up to 22 occasions. Nonmem creates an ETA normal variable for each individual corresponding to each $OMEGA above. Hence, each individual will have 4 ETAs (with 2 set to zero) and an additional 22 occasion ETAS.

You just have to do all the indexing in your head to use the ETAs in the subsequent PK block. In contrast the PoPy script defines the equivalent var f[X] variables as:-

EFFECTS:
    POP: |
        f[TCVL_isv, TVV1_isv] ~ spd_matrix() [
            [ 0.175278 ],
            [ 0.115896, 0.112362 ]
        ]
        f[TVQ_isv] = 0
        f[TVV2_isv] ~ P 0.131759
        f[TVQ2_isv] = 0
        f[TVV3_isv] ~ P 0.177214
        f[TVCL_iov] ~ P 0.0140684

Here the inter-subject variability (isv) variables are just more f[X] variables with convenient labels. Note in contrast to Nonmem the f[TVCL_iov] occasion variance is just defined once here (it being one fixed effect).

In PoPy the mechanism for creating the appropriate number of r[X] variables given the var f[X] is as shown in the ‘ID’ and ‘OCC’ sections:-

EFFECTS:
    ID: |
        r[TVCL_isv, TVV1] ~ mnorm([0,0], f[TCVL_isv, TVV1_isv])
        r[TVQ] ~ norm(0, f[TVQ_isv])
        r[TVV2] ~ norm(0, f[TVV2_isv])
        r[TVQ2] ~ norm(0, f[TVQ2_isv])
        r[TVV3] ~ norm(0, f[TVV3_isv])
    OCC: |
        r[TVCL_iov] ~ norm(0, f[TVCL_iov])

This syntax means that each individual has one r[TVCL_isv, TVV1] and one r[TVQ] instance etc. However the number of instances of r[TVCL_iov] are dependent on the number of values of the ‘OCC’ field for each individual in the data file. i.e. Each individual will have a different sample of r[TVCL_iov] for each occasion. You do not need to know the highest occasion number here (unlike in Nonmem).

In the Nonmem script, the sigma variable definitions are:-

$SIGMA  0.036033  ; variance PROP res error, initial estimate
$SIGMA  0.0164023 ; additional res error, initial estimate

In the PoPy script the equivalent noise f[X] variables are:-

EFFECTS:
    POP: |
        f[PNOISE_var] ~ P 0.036033
        f[ANOISE_var] ~ P 0.0164023

The Nonmem ‘PK’ section is defined as follows:-

$PK
    ; Three-comp model
    IF(NEWIND.NE.2)OTIM1=0
    IF(NEWIND.NE.2)OCOV1=0
    IF(NEWIND.NE.2)OTIM2=0
    IF(NEWIND.NE.2)OCOV2=0
    ;
    STUDY=0
    IF(ID.LT.2000) STUDY=1                     ;Glasgow, Thomson1988
    IF(ID.GE.2000.AND.ID.LT.3000) STUDY=2      ;Uppsala, Nielsen2009
    IF(ID.GE.3000) STUDY=3                     ;Estonia, unpublished
    ;
    WTKG = WT/1000
    ;
    T50  = THETA(7)
    HILL = THETA(8)
    MF   = PMA**HILL/(PMA**HILL+T50**HILL)
    ;
    CREAT2 = CREAT
    IF(CREAT.LT.0) CREAT2 = TCREA       ; when SCr is NA=-99, it is the typical SCr
    ;OF = (CREAT2/TCREA)**(THETA(9))
    ;
    P50  = THETA(10)
    ;PNAF = PNA/(P50+PNA)
    ;
    CRPWR = THETA(9)
    ;IOV code
    BOVC = 0
    IF(OCC.EQ.1)  BOVC = ETA(7)
    IF(OCC.EQ.2)  BOVC = ETA(8)
    IF(OCC.EQ.3)  BOVC = ETA(9)
    IF(OCC.EQ.4)  BOVC = ETA(10)
    IF(OCC.EQ.5)  BOVC = ETA(11)
    IF(OCC.EQ.6)  BOVC = ETA(12)
    IF(OCC.EQ.7)  BOVC = ETA(13)
    IF(OCC.EQ.8)  BOVC = ETA(14)
    IF(OCC.EQ.9)  BOVC = ETA(15)
    IF(OCC.EQ.10) BOVC = ETA(16)
    IF(OCC.EQ.11) BOVC = ETA(17)
    IF(OCC.EQ.12) BOVC = ETA(18)
    IF(OCC.EQ.13) BOVC = ETA(19)
    IF(OCC.EQ.14) BOVC = ETA(20)
    IF(OCC.EQ.15) BOVC = ETA(21)
    IF(OCC.EQ.16) BOVC = ETA(22)
    IF(OCC.EQ.17) BOVC = ETA(23)
    IF(OCC.EQ.18) BOVC = ETA(24)
    IF(OCC.EQ.19) BOVC = ETA(25)
    IF(OCC.EQ.20) BOVC = ETA(26)
    IF(OCC.EQ.21) BOVC = ETA(27)
    IF(OCC.EQ.22) BOVC = ETA(28)
    ;
    TVCL = THETA(1)*MF*(WTKG/70)**(0.632)   ; typical value of CL
    TVV1 = THETA(2)*(WTKG/70)                       ; typical value of V1
    TVQ  = THETA(3)*(WTKG/70)**(0.75)               ; ty. value of intercompartmental CL
    TVV2 = THETA(4)*(WTKG/70)                       ; ty. value of V2
    TVQ2 = THETA(5)*(WTKG/70)**(0.75)               ; ty value of CL3
    TVV3 = THETA(6)*(WTKG/70)                       ; ty value of V3
    ;
    CL   = TVCL*EXP(ETA(1)+BOVC)  ; individual value of CL
    V1   = TVV1*EXP(ETA(2))
    Q    = TVQ*EXP(ETA(3))
    V2   = TVV2*EXP(ETA(4))
    Q2   = TVQ2*EXP(ETA(5))
    V3   = TVV3*EXP(ETA(6))
    ;
    K    = CL/V1
    K12  = Q/V1
    K21  = Q/V2
    K13  = Q2/V1
    K31  = Q2/V3
    ;
    IF(EVID.EQ.1) TM=TIME
    IF(EVID.EQ.1) TAD=0
    IF(EVID.NE.1) TAD=TIME-TM
    ;
    SL1 = 0
    IF(TIME.GT.OTIM1) SL1 = (PNA-OCOV1)/(TIME-OTIM1)
    A_0(4) = PNA
    ;
    SL2 = 0
    IF(TIME.GT.OTIM2) SL2 = (CREAT2-OCOV2)/(TIME-OTIM2)
    A_0(5) = CREAT2

The equivalent in PoPy is the MODEL_PARAMS section:-

MODEL_PARAMS: |

    WTKG = c[WT] / 1000
    MF = c[PMA]**f[HILL]/(c[PMA]**f[HILL] + f[T50]**f[HILL])
    TVCL = f[TVCL] * (WTKG/70)**(0.632) * MF
    TVV1 = f[TVV1] * (WTKG/70)
    TVQ  = f[TVQ]  * (WTKG/70)**(0.75)
    TVV2 = f[TVV2] * (WTKG/70)
    TVQ2 = f[TVQ2] * (WTKG/70)**(0.75)
    TVV3 = f[TVV3] * (WTKG/70)
    BOVC = r[TVCL_iov]
    CL = TVCL * exp(r[TVCL_isv] + BOVC)
    V1 = TVV1 * exp(r[TVV1])
    Q  = TVQ * exp(r[TVQ])
    V2 = TVV2 * exp(r[TVV2])
    Q2 = TVQ2 * exp(r[TVQ2])
    V3 = TVV3 * exp(r[TVV3])
    m[K]   = CL/V1
    m[K12] = Q/V1
    m[K21] = Q/V2
    m[K13] = Q2/V1
    m[K31] = Q2/V3
    m[CRPWR] = f[CRPWR]
    m[P50] = f[P50]
    m[SL1] = 0
    if c[CREAT] < 0:
        m[CREAT2] = c[TCREA]
    else:
        m[CREAT2] = c[CREAT]
    m[SL2] = 0
    m[V1] = V1
    m[PNOISE_var] = f[PNOISE_var]
    m[ANOISE_var] = f[ANOISE_var]

The PoPy version is shorter, mainly because it does not have to write an if statement for all values of the OCC variable i.e. ‘IF(OCC.EQ.1) BOVC = ETA(7)’ etc.

Note that the Nonmem $PK section contains the lines:-

A_0(4) = PNA
A_0(5) = CREAT2

These are necessary to provide initial values for the $DES section. The equivalent in PoPy is the STATES section:-

STATES: |
    s[COVCMT1] = c[PNA]
    s[COVCMT2] = m[CREAT2]

which provides s[X] values at t=0.0 for the PoPy DERIVATIVES section. Note by default in both PoPy s[X] =0.0, unless the STATES section says otherwise. A similar convention is used in Nonmem.

The Nonmem $DES section is as follows:-

$DES
    TCOV1 = A(4)
    TCOV2 = A(5)
    PNAF = TCOV1/(P50+TCOV1)
    OF = (TCOV2/TCREA)**CRPWR

    DADT(1) = A(3)*K31+A(2)*K21-A(1)*(K*PNAF*OF+K12+K13)
    DADT(2) = A(1)*K12-A(2)*K21
    DADT(3) = A(1)*K13-A(3)*K31
    DADT(4)= SL1
    DADT(5)= SL2

This is very similar to the PoPy DERIVATIVES section:-

DERIVATIVES: |
    TCOV1 = s[COVCMT1]
    TCOV2 = s[COVCMT2]
    PNAF = TCOV1 / (m[P50]+TCOV1)
    OF = (TCOV2/c[TCREA])**m[CRPWR]
    dose[DOSE_infrate] = @inf_rate{ amt: c[AMT], rate: c[RATE], lag: 0 }
    d[CENTRAL] = dose[DOSE_infrate] + m[K31]*s[PERIPH2] + m[K21]*s[PERIPH1] - (m[K]*PNAF*OF+m[K12]+m[K13])*s[CENTRAL]
    d[PERIPH1] = - m[K21]*s[PERIPH1] + m[K12] *s[CENTRAL]
    d[PERIPH2] = - m[K31]*s[PERIPH2] + m[K13] *s[CENTRAL]
    d[COVCMT1] =   m[SL1]
    d[COVCMT2] =   m[SL2]

except that in PoPy the dose[DOSE_infrate] is defined and explicitly placed in the s[CENTRAL] compartment. Nonmem puts all boluses/infusions in compartment one by convention, unless a CMT field is defined in the data set (it isn’t in this case). In Nonmem you have to examine the data file to determine the type of dose, however in the PoPy syntax above it’s pretty clear it’s an infusion defined by the c[AMT] and c[RATE] data entries in each dosing data row.

The ordinary differential equation solver parameters for Nonmem are set here:-

$SUBROUTINE ADVAN6 TOL=6

Here ‘ADVAN6’ is the Nonmem ode solver for non-stiff systems. With a relative tolerance of 1e-6. Note this solver is different from the CPPLSODA method employed by PoPy which is similar to selecting ‘ADVAN13’ in Nonmem. In this case the ordinary differential equation solver method makes little difference, as both PoPy and Nonmem return similar objective values for this model.

In PoPy the ordinary differential equation solver parameters are defined as follows:-

ODE_SOLVER:
    CPPLSODA:
        # Absolute tolerance of ode solver.
        atol: 1e-12

        # Relative tolerance of ode solver.
        rtol: 1e-12

        # Maximum number of steps allowed in ode solver.
        max_nsteps: 10000000

The $ERROR block in the Nonmem script defines a proportional + additional noise model:-

$ERROR
    IPRED  = A(1)/V1
    Y      = IPRED*(1+EPS(1)) + EPS(2)

The PREDICTIONS section in the PoPy script defines the same error model:-

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

The Nonmem script defines the FOCE fitting method (METHOD=1) here:-

$ESTIMATION METHOD=1 INTER MAXEVAL=0 PRINT=1

The equivalent in the PoPy script is here:-

FIT_METHODS: [ND: {max_n_main_iterations: 0}]

Note both Nonmem and PoPy only optimize the r[X] in the fitting run as the number of main iterations is set to zero.