# DDMoRe0238 Conversion Example

## 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

Here we describe the manual conversion of the Nonmem script file, `c:\PoPy\validation\ddmore0238\Executable_run35b_ddm2.ctl` to create the PoPy Fit Script `c:\PoPy\validation\ddmore0238\popy\fit_script.pyml`.

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

### Data Source

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: Simulated_simdataDDM.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 57.

Table 57 Main data fields for Nonmem (first eight rows)
ID | TIME | RATE | EVID | AMT | WT | CREAT | DV | OCC
1001 | 0 | 72 | 1 | 6 | 2120 | 78 | 0 | 1
1001 | 11.5 | 0 | 0 | 0 | 2120 | 78 | 1.109 | 1
1001 | 12 | 72 | 1 | 6 | 2120 | 78 | 0 | 2
1001 | 12.5 | 0 | 0 | 0 | 2120 | 78 | 7.0181 | 2
1001 | 24 | 48 | 1 | 4 | 2120 | 78 | 0 | 2
1001 | 36 | 48 | 1 | 4 | 2120 | 78 | 0 | 2
1001 | 48 | 48 | 1 | 4 | 2120 | 78 | 0 | 3
1001 | 58.5 | 0 | 0 | 0 | 2120 | 78 | 2.2136 | 3

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 indicating that the dose is an infusion (RATE>0).
        if c[AMT] > 0 and c[RATE] > 0:
            c[TYPE] += ":DOSE_infrate"

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

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

Table 58 Extra data fields for PoPy (first eight rows)
DV_CENTRAL | DV_CENTRAL_FLAG | TYPE
0 | 0 | dose:DOSE_infrate
1.109 | 1 | obs
0 | 0 | dose:DOSE_infrate
7.0181 | 1 | obs
0 | 0 | dose:DOSE_infrate
0 | 0 | dose:DOSE_infrate
0 | 0 | dose:DOSE_infrate
2.2136 | 1 | obs

### Fixed/Random Effects Specification

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

### PK Model Parameters Specification

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.

### Initial State Specification

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.

### ODEs Specification

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

### Likelihood Error Specification

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)

### Parameter Fitting Methods Specification

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.