# DDMoRe0061 Conversion Example

## Overview

A real world model, based on a gentamicin PK study [Harling2015]. It has a combined bolus + infusion dose and the original PK model uses ADVAN3 TRANS4 in Nonmem. There are some covariate effects on the clearance and volume of distribution of the central compartment parameters. There are 210 individuals and 1949 data rows. See DDMoRe: 0061

Here we describe the manual conversion of the Nonmem script file `c:\PoPy\validation\ddmore0061\Executable_gentamicin_pk.mod` to create the PoPy Fit Script `c:\PoPy\validation\ddmore0061\popy\fit_script.pyml`.

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

Table 55 Main data fields for Nonmem (first six rows)
ID | TIME | AMT | RATE | DUR | DV | MDV | EVID
1 | 0 | 150 | 0 | 0 | 12.376 | 1 | 1
1 | 9 | 150 | 0 | 0 | 13.156 | 1 | 1
1 | 14.25 | 0 | 0 | 0 | 2.2347 | 0 | 0
1 | 16.92 | 150 | 0 | 0 | 13.599 | 1 | 1
1 | 25 | 150 | 250 | 0.5 | 1.322 | 1 | 1
1 | 33 | 0 | 0 | 0 | 1.459 | 0 | 0

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.

Table 56 Extra data fields for PoPy (first six rows)
DRUG_CONC | DRUG_CONC_FLAG | TYPE
0 | 0 | dose:DOSE_bolus
0 | 0 | dose:DOSE_bolus
2.2347 | 1 | obs
0 | 0 | dose:DOSE_bolus
0 | 0 | dose:DOSE_infrate
1.459 | 1 | obs

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.