# DDMoRe0061 Conversion Example

## Overview

A real world model, based on a gentamicin PK study [Harling2015]. Combined bolus + infusion PK model. Uses ADVAN3 TRANS4 in Nonmem. Some covariate effects on PK parameters. 210 individuals. 1949 total data rows. See DDMoRe: 0061

## Data Conversion

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

    c:\PoPy\validation\ddmore0061\Simulated_gentamicin_pk.csv

To this data set in PoPy format:-

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

Using the N2PDat Script located here:-

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

The original Nonmem data set looks something like Table 48.

Table 48 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

The critical entries for the N2PDat Script are as follows:-

    INPUT_NONMEM_FIELDS:
         obs_cmt_numbers: [1]
         dose_cmt_numbers: [1]

These two entries specify that the observations and doses all occur in compartment one. Technically the Nonmem data set above has no CMT field. However CMT defaults to one in Nonmem if it is NOT provided.

The dataset output by N2PDat Script is the same as the input, but with the 3 extra columns shown in Table 49.

Table 49 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 the Nonmem ‘MDV’ flag. The new ‘TYPE’ field is formatted to process either bolus or infusion doses in PoPy.

## Script Conversion

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

    c:\PoPy\validation\ddmore0061\Executable_gentamicin_pk.mod

To create the PoPy Fit Script here:-

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

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: popy_data.csv}

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

    LEVEL_PARAMS:
        GLOBAL:
            params: |
                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 to be clumsily 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:-

    LEVEL_PARAMS:
        INDIV:
            params: |
                r[CL_isv, Q_isv] ~ mnorm([0,0], f[CL_isv, Q_isv])
            split_field: ID
            split_dict: "*"

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 sigmas are fixed effects and also defined in LEVEL_PARAMS.

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 the PoPy 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.

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: {SCIPY_ODEINT: {}}

Which tells PoPy to use the scipy ‘odeint’ solver, which is based on LSODA (like Nonmem Advan13).

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. It is not clear (to the PoPy developers) 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.

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 JOE fitting method (roughly equivalent to FOCE) is:-

    FIT_METHODS: [JOE: {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.