Modeling

While the pharmpy.model.Model class can be directly manipulated with low level operations the modeling module offers higher level operations and transformations for building a model. These transformations are also available via the Pharmpy command line interface. To read more about these functions such as how the initial estimates of parameters are chosen, see their respective API documentation.

Basic functions

Reading and writing models

Pharmpy can read in NONMEM models using the pharmpy.modeling.read_model() function:

PythonR

from pharmpy.modeling import *
model = read_model(path / 'pheno.mod')

To inspect the read in model code, the pharmpy.modeling.print_model_code() function is available:

PythonR

print_model_code(model)

If the model code is in a string variable it can be read in directly using pharmpy.modeling.read_model_from_string().

PythonR

code = '$PROBLEM base model\n$INPUT ID DV TIME\n$DATA file.csv IGNORE=@\n$PRED Y = THETA(1) + ETA(1) + ERR(1)\n$THETA 0.1\n$OMEGA 0.01\n$SIGMA 1\n$ESTIMATION METHOD=1'
model = read_model_from_string(code)

Finally, any Pharmpy model can be written to a file with pharmpy.modeling.write_model():

PythonR

write_model(model, 'mymodel.mod')

Loading example models

Pharmpy has example models with example datasets that can be accessed using pharmpy.modeling.load_example_model():

PythonR

model = load_example_model('pheno')
print_model_code(model)

$PROBLEM PHENOBARB SIMPLE MODEL
$DATA pheno.dta IGNORE=@
$INPUT ID TIME AMT WGT APGR DV FA1 FA2
$SUBROUTINE ADVAN1 TRANS2
$ABBREV REPLACE ETA_CL=ETA(1)
$ABBREV REPLACE ETA_VC=ETA(2)

$PK
TVCL = THETA(1)*WGT
TVV = THETA(2)*WGT
IF(APGR.LT.5) TVV = TVV*(1 + THETA(3))
CL = TVCL*EXP(ETA_CL)
VC = TVV*EXP(ETA_VC)
V = VC
S1 = VC

$ERROR
Y = F + F*EPS(1)

$THETA  (0,0.00469307) ; POP_CL
$THETA  (0,1.00916) ; POP_VC
$THETA  (-.99,.1) ; COVAPGR

$OMEGA  0.0309626 ; IIV_CL
$OMEGA  0.031128 ; IIV_VC

$SIGMA  0.0130865  ; SIGMA

$ESTIMATION METHOD=1 INTERACTION MAXEVALS=99999
$COVARIANCE UNCONDITIONAL PRINT=E
$TABLE ID TIME DV CIPREDI PRED RES CWRES NOAPPEND NOPRINT ONEHEADER FILE=pheno.tab

Converting models

Pharmpy supports the estimation software NONMEM, nlmixr2 and rxode2, and a Pharmpy model can be converted into those formats using pharmpy.modeling.convert_model():

PythonR

model_nlmixr = convert_model(model, 'nlmixr')

Then we can inspect the new model code:

PythonR

print_model_code(model_nlmixr)

pheno <- function() {
    ini({
        # --- THETAS ---
        POP_CL <- c(0.0, 0.00469307, Inf)
        POP_VC <- c(0.0, 1.00916, Inf)
        COVAPGR <- c(-0.99, 0.1, Inf)
        
        # --- ETAS ---
        ETA_CL ~ 0.0309626
        ETA_VC ~ 0.031128
        
        # --- EPSILONS ---
        SIGMA <- c(0.0, 0.0130865, Inf)
    })
    model({
        TVCL <- POP_CL*WGT
        TVV <- POP_VC*WGT
        if (APGR < 5) {
            TVV <- TVV*(COVAPGR + 1)
        } else {
            TVV <- TVV
        }
        CL <- TVCL*exp(ETA_CL)
        VC <- TVV*exp(ETA_VC)
        V <- VC
        S1 <- VC
        
        # --- DIFF EQUATIONS ---
        d/dt(A_CENTRAL) = -CL*A_CENTRAL/V
        
        F <- A_CENTRAL/S1
        Y <- F
        add_error <- 0
        prop_error <- SIGMA
        Y ~ add(add_error) + prop(prop_error)
    })
}

fit <- nlmixr2(pheno, dataset, est = "focei", control=foceiControl(maxOuterIterations=99999))

Create basic models

As a starting point for this user guide, we will create a basic PK model using the function pharmpy.modeling.create_basic_pk_model():

PythonR

dataset_path = path / 'pheno.dta'
model_start = create_basic_pk_model(administration='oral',
                                    dataset_path=dataset_path,
                                    cl_init=0.01,
                                    vc_init=1.0,
                                    mat_init=0.1)

We can examine the model statements of the model:

PythonR

model_start.statements

\begin{align*}MAT &= POP_{MAT} \cdot e^{\eta_{MAT}}\\CL &= POP_{CL} \cdot e^{\eta_{CL}}\\VC &= POP_{VC} \cdot e^{\eta_{VC}}\\\end{align*}

Bolus(AMT, admid=1) → DEPOT 
┌─────┐        ┌───────┐        
│DEPOT│──1/MAT→│CENTRAL│──CL/VC→
└─────┘        └───────┘

\begin{align*}IPRED &= \frac{A_{CENTRAL}{\left(t \right)}}{VC}\\IPREDADJ &= \begin{cases} 2.225 \cdot 10^{-16} & \text{for}\: IPRED &= 0 \\IPRED & \text{otherwise} \end{cases}\\Y &= IPRED + IPREDADJ \cdot \epsilon_{p}\\\end{align*}

We can see that the model is a one compartment model with first order absorption and elimination and no absorption delay. Examining the random variables:

PythonR

model_start.random_variables

\[\begin{split}\begin{align*} \left[\begin{matrix}\eta_{CL}\\\eta_{VC}\end{matrix}\right] & \sim \mathcal{N} \left(\left[\begin{matrix}0\\0\end{matrix}\right],\left[\begin{matrix}IIV_{CL} & IIV_{CL IIV VC}\\IIV_{CL IIV VC} & IIV_{VC}\end{matrix}\right]\right) \\ \epsilon_{p} & \sim \mathcal{N} \left(0,\sigma\right) \\ \eta_{MAT} & \sim \mathcal{N} \left(0,IIV_{MAT}\right)\end{align*}\end{split}\]

Next we can convert the start model from a generic Pharmpy model to a NONMEM model:

PythonR

model_start = convert_model(model_start, 'nonmem')

We can then examine the NONMEM model code:

PythonR

print_model_code(model_start)

$PROBLEM Start model
$INPUT ID TIME AMT WGT APGR DV FA1 FA2
$DATA file.csv IGNORE=@
$SUBROUTINES ADVAN2 TRANS2
$ABBR REPLACE ETA_CL=ETA(1)
$ABBR REPLACE ETA_VC=ETA(2)
$ABBR REPLACE ETA_MAT=ETA(3)
$PK
MAT = THETA(3)*EXP(ETA_MAT)
CL = THETA(1)*EXP(ETA_CL)
VC = THETA(2)*EXP(ETA_VC)
KA = 1/MAT
V = VC
$ERROR
IPRED = A(2)/VC
IF (IPRED.EQ.0) THEN
    IPREDADJ = 2.22500000000000E-16
ELSE
    IPREDADJ = IPRED
END IF
Y = IPRED + EPS(1)*IPREDADJ
$THETA  (0,0.01) ; POP_CL
$THETA  (0,1.0) ; POP_VC
$THETA  (0,0.1) ; POP_MAT
$OMEGA BLOCK(2)
0.1	; IIV_CL
0.01	; IIV_CL_IIV_VC
0.1	; IIV_VC
$OMEGA  0.1 ; IIV_MAT
$SIGMA  0.09 ; sigma
$TABLE ID TIME DV CIPREDI PRED CWRES FILE=mytab ONEHEADER NOAPPEND NOPRINT
$ESTIMATION METHOD=COND INTER MAXEVAL=99999

Modeling transformations

In Pharmpy there are many different modeling functions that modify the model object. In Pharmpy, the model object and all its attributes are immutable, meaning that the modeling functions always return a copy of the model object.

Note

To see more information on how initial estimates are chosen etc., please check the API reference.

Structural model

There are many functions to change or examine the structural model of a PK dataset. Using the base model from above, we’ll go through how to change different aspects of the structural model step by step.

Absorption rate

As an example, we’ll set the absorption of the start model to zero order absorption:

PythonR

run1 = set_zero_order_absorption(model_start)
run1.statements.ode_system

Infusion(AMT, admid=1, duration=D1) → CENTRAL 
┌───────┐       
│CENTRAL│──CL/V→
└───────┘

And examine the updated NONMEM code:

PythonR

print_model_code(run1)

$PROBLEM Start model
$INPUT ID TIME AMT WGT APGR DV FA1 FA2 RATE
$DATA DUMMYPATH IGNORE=@
$SUBROUTINES ADVAN1 TRANS2
$ABBR REPLACE ETA_CL=ETA(1)
$ABBR REPLACE ETA_VC=ETA(2)
$ABBR REPLACE ETA_MAT=ETA(3)
$PK
MAT = THETA(3)*EXP(ETA_MAT)
CL = THETA(1)*EXP(ETA_CL)
VC = THETA(2)*EXP(ETA_VC)
V = VC
D1 = 2*MAT
$ERROR
IPRED = A(1)/VC
IF (IPRED.EQ.0) THEN
    IPREDADJ = 2.22500000000000E-16
ELSE
    IPREDADJ = IPRED
END IF
Y = IPRED + EPS(1)*IPREDADJ
$THETA  (0,0.01) ; POP_CL
$THETA  (0,1.0) ; POP_VC
$THETA  (0,0.1) ; POP_MAT
$OMEGA BLOCK(2)
0.1	; IIV_CL
0.01	; IIV_CL_IIV_VC
0.1	; IIV_VC
$OMEGA  0.1 ; IIV_MAT
$SIGMA  0.09 ; sigma
$TABLE ID TIME DV CIPREDI PRED CWRES FILE=mytab ONEHEADER NOAPPEND NOPRINT
$ESTIMATION METHOD=COND INTER MAXEVAL=99999

Note that the ADVAN has been updated.

List of functions to change absorption rate:

• pharmpy.modeling.set_bolus_absorption()

• pharmpy.modeling.set_first_order_absorption()

• pharmpy.modeling.set_seq_zo_fo_absorption()

• pharmpy.modeling.set_zero_order_absorption()

Absorption delay

Next, we will add absorption delay.

PythonR

run2 = add_lag_time(run1)
run2.statements.ode_system

Infusion(AMT, admid=1, duration=D1) → CENTRAL 
┌──────────┐       
│ CENTRAL  │       
│tlag=ALAG1│──CL/V→
└──────────┘

And examine the model code:

PythonR

print_model_code(run2)

$PROBLEM Start model
$INPUT ID TIME AMT WGT APGR DV FA1 FA2 RATE
$DATA DUMMYPATH IGNORE=@
$SUBROUTINES ADVAN1 TRANS2
$ABBR REPLACE ETA_CL=ETA(1)
$ABBR REPLACE ETA_VC=ETA(2)
$ABBR REPLACE ETA_MAT=ETA(3)
$PK
MDT = THETA(4)
MAT = THETA(3)*EXP(ETA_MAT)
CL = THETA(1)*EXP(ETA_CL)
VC = THETA(2)*EXP(ETA_VC)
V = VC
D1 = 2*MAT
ALAG1 = MDT
$ERROR
IPRED = A(1)/VC
IF (IPRED.EQ.0) THEN
    IPREDADJ = 2.22500000000000E-16
ELSE
    IPREDADJ = IPRED
END IF
Y = IPRED + EPS(1)*IPREDADJ
$THETA  (0,0.01) ; POP_CL
$THETA  (0,1.0) ; POP_VC
$THETA  (0,0.1) ; POP_MAT
$THETA  (0,0.5) ; POP_MDT
$OMEGA BLOCK(2)
0.1	; IIV_CL
0.01	; IIV_CL_IIV_VC
0.1	; IIV_VC
$OMEGA  0.1 ; IIV_MAT
$SIGMA  0.09 ; sigma
$TABLE ID TIME DV CIPREDI PRED CWRES FILE=mytab ONEHEADER NOAPPEND NOPRINT
$ESTIMATION METHOD=COND INTER MAXEVAL=99999

List of functions to change elimination:

• pharmpy.modeling.add_lag_time()

• pharmpy.modeling.remove_lag_time()

• pharmpy.modeling.set_transit_compartments()

Distribution

It is possible to change the number of peripheral compartments. Let us add one peripheral compartment

PythonR

run3 = add_peripheral_compartment(run2)
run3.statements.ode_system

Infusion(AMT, admid=1, duration=D1) → CENTRAL 
┌───────────┐        
│PERIPHERAL1│        
└───────────┘        
  ↑       │          
Q/V1     Q/V2        
  │       ↓          
┌───────────┐        
│  CENTRAL  │        
│tlag=ALAG1 │──CL/V1→
└───────────┘

And examine the model code:

PythonR

print_model_code(run3)

$PROBLEM Start model
$INPUT ID TIME AMT WGT APGR DV FA1 FA2 RATE
$DATA DUMMYPATH IGNORE=@
$SUBROUTINES ADVAN3 TRANS4
$ABBR REPLACE ETA_CL=ETA(1)
$ABBR REPLACE ETA_VC=ETA(2)
$ABBR REPLACE ETA_MAT=ETA(3)
$PK
VP1 = THETA(6)
QP1 = THETA(5)
MDT = THETA(4)
MAT = THETA(3)*EXP(ETA_MAT)
CL = THETA(1)*EXP(ETA_CL)
VC = THETA(2)*EXP(ETA_VC)
V1 = VC
D1 = 2*MAT
ALAG1 = MDT
Q = QP1
V2 = VP1
$ERROR
IPRED = A(1)/VC
IF (IPRED.EQ.0) THEN
    IPREDADJ = 2.22500000000000E-16
ELSE
    IPREDADJ = IPRED
END IF
Y = IPRED + EPS(1)*IPREDADJ
$THETA  (0,0.01) ; POP_CL
$THETA  (0,1.0) ; POP_VC
$THETA  (0,0.1) ; POP_MAT
$THETA  (0,0.5) ; POP_MDT
$THETA  (0,0.01) ; POP_QP1
$THETA  (0,0.05) ; POP_VP1
$OMEGA BLOCK(2)
0.1	; IIV_CL
0.01	; IIV_CL_IIV_VC
0.1	; IIV_VC
$OMEGA  0.1 ; IIV_MAT
$SIGMA  0.09 ; sigma
$TABLE ID TIME DV CIPREDI PRED CWRES FILE=mytab ONEHEADER NOAPPEND NOPRINT
$ESTIMATION METHOD=COND INTER MAXEVAL=99999

List of functions to change distribution:

• pharmpy.modeling.add_peripheral_compartment()

• pharmpy.modeling.remove_peripheral_compartment()

• pharmpy.modeling.set_peripheral_compartments()

Elimination

Now we will change to non-linear elimination.

PythonR

run4 = set_michaelis_menten_elimination(run3)
run4.statements.ode_system

Infusion(AMT, admid=1, duration=D1) → CENTRAL 
┌───────────┐                                      
│PERIPHERAL1│                                      
└───────────┘                                      
  ↑       │                                        
Q/V1     Q/V2                                      
  │       ↓                                        
┌───────────┐                                      
│  CENTRAL  │                                      
│tlag=ALAG1 │──CLMM*KM/(V1*(KM + A_CENTRAL(t)/V1))→
└───────────┘

And examine the model code:

PythonR

print_model_code(run4)

$PROBLEM Start model
$INPUT ID TIME AMT WGT APGR DV FA1 FA2 RATE
$DATA DUMMYPATH IGNORE=@
$SUBROUTINES ADVAN13 TOL=9
$MODEL COMPARTMENT=(CENTRAL DEFDOSE) COMPARTMENT=(PERIPHERAL1)
$ABBR REPLACE ETA_CL=ETA(1)
$ABBR REPLACE ETA_VC=ETA(2)
$ABBR REPLACE ETA_MAT=ETA(3)
$PK
KM = THETA(7)
VP1 = THETA(6)
QP1 = THETA(5)
MDT = THETA(4)
MAT = THETA(3)*EXP(ETA_MAT)
CLMM = THETA(1)*EXP(ETA_CL)
VC = THETA(2)*EXP(ETA_VC)
V1 = VC
D1 = 2*MAT
ALAG1 = MDT
Q = QP1
V2 = VP1
$DES
DADT(1) = A(1)*(-CLMM*KM/(V1*(A(1)/V1 + KM)) - Q/V1) + A(2)*Q/V2
DADT(2) = A(1)*Q/V1 - A(2)*Q/V2
$ERROR
IPRED = A(1)/VC
IF (IPRED.EQ.0) THEN
    IPREDADJ = 2.22500000000000E-16
ELSE
    IPREDADJ = IPRED
END IF
Y = IPRED + EPS(1)*IPREDADJ
$THETA  (0,0.01) ; POP_CLMM
$THETA  (0,1.0) ; POP_VC
$THETA  (0,0.1) ; POP_MAT
$THETA  (0,0.5) ; POP_MDT
$THETA  (0,0.01) ; POP_QP1
$THETA  (0,0.05) ; POP_VP1
$THETA  (0,33.95,101.85000000000001) ; POP_KM
$OMEGA BLOCK(2)
0.1	; IIV_CLMM
0.01	; IIV_CLMM_IIV_VC
0.1	; IIV_VC
$OMEGA  0.1 ; IIV_MAT
$SIGMA  0.09 ; sigma
$TABLE ID TIME DV CIPREDI PRED CWRES FILE=mytab ONEHEADER NOAPPEND NOPRINT
$ESTIMATION METHOD=COND INTER MAXEVAL=99999

List of functions to change elimination:

• pharmpy.modeling.set_first_order_elimination()

• pharmpy.modeling.set_michaelis_menten_elimination()

• pharmpy.modeling.set_mixed_mm_fo_elimination()

• pharmpy.modeling.set_zero_order_elimination()

Parameter variability model

Pharmpy has multiple functions to change the parameter variability model. Using the base model from above, we’ll go through different aspects of changing the parameter variability model.

Adding and removing parameter variability

It is possible to add and remove inter-individual variability (IIV) and inter-occasion variability (IOV) using pharmpy.modeling.add_iiv() and pharmpy.modeling.add_iov(). Since the start model from above has IIV on all its parameters, we will start by removing an IIV using the pharmpy.modeling.remove_iiv() function.

PythonR

run1 = remove_iiv(model_start, 'MAT')
run1.random_variables.iiv

\[\begin{split}\begin{align*} \left[\begin{matrix}\eta_{CL}\\\eta_{VC}\end{matrix}\right] & \sim \mathcal{N} \left(\left[\begin{matrix}0\\0\end{matrix}\right],\left[\begin{matrix}IIV_{CL} & IIV_{CL IIV VC}\\IIV_{CL IIV VC} & IIV_{VC}\end{matrix}\right]\right)\end{align*}\end{split}\]

Next, we add an IIV to the same parameter:

PythonR

run2 = add_iiv(run1, 'MAT', 'exp', operation='*')
run2.random_variables.iiv

\[\begin{split}\begin{align*} \left[\begin{matrix}\eta_{CL}\\\eta_{VC}\end{matrix}\right] & \sim \mathcal{N} \left(\left[\begin{matrix}0\\0\end{matrix}\right],\left[\begin{matrix}IIV_{CL} & IIV_{CL IIV VC}\\IIV_{CL IIV VC} & IIV_{VC}\end{matrix}\right]\right) \\ \eta_{MAT} & \sim \mathcal{N} \left(0,IIV_{MAT}\right)\end{align*}\end{split}\]

And examine the model code:

PythonR

print_model_code(run2)

$PROBLEM Start model
$INPUT ID TIME AMT WGT APGR DV FA1 FA2
$DATA file.csv IGNORE=@
$SUBROUTINES ADVAN2 TRANS2
$ABBR REPLACE ETA_CL=ETA(1)
$ABBR REPLACE ETA_VC=ETA(2)
$ABBR REPLACE ETA_MAT=ETA(3)
$PK
MAT = THETA(3)*EXP(ETA_MAT)
CL = THETA(1)*EXP(ETA_CL)
VC = THETA(2)*EXP(ETA_VC)
KA = 1/MAT
V = VC
$ERROR
IPRED = A(2)/VC
IF (IPRED.EQ.0) THEN
    IPREDADJ = 2.22500000000000E-16
ELSE
    IPREDADJ = IPRED
END IF
Y = IPRED + EPS(1)*IPREDADJ
$THETA  (0,0.01) ; POP_CL
$THETA  (0,1.0) ; POP_VC
$THETA  (0,0.1) ; POP_MAT
$OMEGA BLOCK(2)
0.1	; IIV_CL
0.01	; IIV_CL_IIV_VC
0.1	; IIV_VC
$OMEGA  0.09 ; IIV_MAT
$SIGMA  0.09 ; sigma
$TABLE ID TIME DV CIPREDI PRED CWRES FILE=mytab ONEHEADER NOAPPEND NOPRINT
$ESTIMATION METHOD=COND INTER MAXEVAL=99999

List of functions to add and remove parameter variability:

• pharmpy.modeling.add_iiv()

• pharmpy.modeling.add_iov()

• pharmpy.modeling.add_pk_iiv()

• pharmpy.modeling.remove_iiv()

• pharmpy.modeling.remove_iov()

Adding and removing covariance

As the next example, we will create a joint distribution using pharmpy.modeling.create_joint_distribution() where the eta on MAT is included:

PythonR

run3 = create_joint_distribution(run2)
run3.random_variables.iiv

\[\begin{split}\begin{align*} \left[\begin{matrix}\eta_{CL}\\\eta_{VC}\\\eta_{MAT}\end{matrix}\right] & \sim \mathcal{N} \left(\left[\begin{matrix}0\\0\\0\end{matrix}\right],\left[\begin{matrix}IIV_{CL} & IIV_{CL IIV VC} & IIV_{CL IIV MAT}\\IIV_{CL IIV VC} & IIV_{VC} & IIV_{VC IIV MAT}\\IIV_{CL IIV MAT} & IIV_{VC IIV MAT} & IIV_{MAT}\end{matrix}\right]\right)\end{align*}\end{split}\]

And examine the model code:

PythonR

print_model_code(run3)

$PROBLEM Start model
$INPUT ID TIME AMT WGT APGR DV FA1 FA2
$DATA file.csv IGNORE=@
$SUBROUTINES ADVAN2 TRANS2
$ABBR REPLACE ETA_CL=ETA(1)
$ABBR REPLACE ETA_VC=ETA(2)
$ABBR REPLACE ETA_MAT=ETA(3)
$PK
MAT = THETA(3)*EXP(ETA_MAT)
CL = THETA(1)*EXP(ETA_CL)
VC = THETA(2)*EXP(ETA_VC)
KA = 1/MAT
V = VC
$ERROR
IPRED = A(2)/VC
IF (IPRED.EQ.0) THEN
    IPREDADJ = 2.22500000000000E-16
ELSE
    IPREDADJ = IPRED
END IF
Y = IPRED + EPS(1)*IPREDADJ
$THETA  (0,0.01) ; POP_CL
$THETA  (0,1.0) ; POP_VC
$THETA  (0,0.1) ; POP_MAT
$OMEGA BLOCK(3)
0.1	; IIV_CL
0.01	; IIV_CL_IIV_VC
0.1	; IIV_VC
0.0094868	; IIV_CL_IIV_MAT
0.0094868	; IIV_VC_IIV_MAT
0.09	; IIV_MAT
$SIGMA  0.09 ; sigma
$TABLE ID TIME DV CIPREDI PRED CWRES FILE=mytab ONEHEADER NOAPPEND NOPRINT
$ESTIMATION METHOD=COND INTER MAXEVAL=99999

List of functions to change covariance structure:

• pharmpy.modeling.create_joint_distribution()

• pharmpy.modeling.split_joint_distribution()

Eta transformations

It is also possible to transform the etas using the following functions:

• pharmpy.modeling.transform_etas_boxcox()

• pharmpy.modeling.transform_etas_john_draper()

• pharmpy.modeling.transform_etas_tdist()

Covariates and allometry

Covariate effects may be applied to a model using the pharmpy.modeling.add_covariate_effect().

PythonR

run1 = add_covariate_effect(model_start, 'CL', 'WGT', 'pow', operation='*')

Here, CL indicates the name of the parameter onto which you want to apply the effect, WGT is the name of the covariate, and pow (power function) is the effect you want to apply. The effect can be either added or multiplied to the parameter, denoted by ‘*’ or ‘+’ (multiplied is default). We can examine the model code:

PythonR

print_model_code(run1)

$PROBLEM Start model
$INPUT ID TIME AMT WGT APGR DV FA1 FA2
$DATA file.csv IGNORE=@
$SUBROUTINES ADVAN2 TRANS2
$ABBR REPLACE ETA_CL=ETA(1)
$ABBR REPLACE ETA_VC=ETA(2)
$ABBR REPLACE ETA_MAT=ETA(3)
$PK
WGT_MEDIAN = 1.30000000000000
MAT = THETA(3)*EXP(ETA_MAT)
CL = THETA(1)*EXP(ETA_CL)
CLWGT = (WGT/WGT_MEDIAN)**THETA(4)
CL = CL*CLWGT
VC = THETA(2)*EXP(ETA_VC)
KA = 1/MAT
V = VC
$ERROR
IPRED = A(2)/VC
IF (IPRED.EQ.0) THEN
    IPREDADJ = 2.22500000000000E-16
ELSE
    IPREDADJ = IPRED
END IF
Y = IPRED + EPS(1)*IPREDADJ
$THETA  (0,0.01) ; POP_CL
$THETA  (0,1.0) ; POP_VC
$THETA  (0,0.1) ; POP_MAT
$THETA  (-100,0.001,100000) ; POP_CLWGT
$OMEGA BLOCK(2)
0.1	; IIV_CL
0.01	; IIV_CL_IIV_VC
0.1	; IIV_VC
$OMEGA  0.1 ; IIV_MAT
$SIGMA  0.09 ; sigma
$TABLE ID TIME DV CIPREDI PRED CWRES FILE=mytab ONEHEADER NOAPPEND NOPRINT
$ESTIMATION METHOD=COND INTER MAXEVAL=99999

Pharmpy also supports user formatted covariate effects.

PythonR

user_effect = '((cov/std) - median) * theta'
run2 = add_covariate_effect(model_start, 'CL', 'WGT', user_effect, operation='*')

The covariate is denoted as cov, the theta as theta (or, if multiple thetas: theta1, theta2 etc.), and the mean, median, and standard deviation as mean, median, and std respectively. This is in order for the names to be substituted with the correct symbols.

PythonR

print_model_code(run2)

$PROBLEM Start model
$INPUT ID TIME AMT WGT APGR DV FA1 FA2
$DATA file.csv IGNORE=@
$SUBROUTINES ADVAN2 TRANS2
$ABBR REPLACE ETA_CL=ETA(1)
$ABBR REPLACE ETA_VC=ETA(2)
$ABBR REPLACE ETA_MAT=ETA(3)
$PK
WGT_MEDIAN = 1.30000000000000
WGT_STD = 0.704564727537012
MAT = THETA(3)*EXP(ETA_MAT)
CL = THETA(1)*EXP(ETA_CL)
CLWGT = THETA(4)*(WGT/WGT_STD - WGT_MEDIAN)
CL = CL*CLWGT
VC = THETA(2)*EXP(ETA_VC)
KA = 1/MAT
V = VC
$ERROR
IPRED = A(2)/VC
IF (IPRED.EQ.0) THEN
    IPREDADJ = 2.22500000000000E-16
ELSE
    IPREDADJ = IPRED
END IF
Y = IPRED + EPS(1)*IPREDADJ
$THETA  (0,0.01) ; POP_CL
$THETA  (0,1.0) ; POP_VC
$THETA  (0,0.1) ; POP_MAT
$THETA  (-100000,0.001,100000) ; POP_CLWGT
$OMEGA BLOCK(2)
0.1	; IIV_CL
0.01	; IIV_CL_IIV_VC
0.1	; IIV_VC
$OMEGA  0.1 ; IIV_MAT
$SIGMA  0.09 ; sigma
$TABLE ID TIME DV CIPREDI PRED CWRES FILE=mytab ONEHEADER NOAPPEND NOPRINT
$ESTIMATION METHOD=COND INTER MAXEVAL=99999

List of functions for covariates and allometry:

• pharmpy.modeling.add_allometry()

• pharmpy.modeling.add_covariate_effect()

• pharmpy.modeling.remove_covariate_effect()

Population parameters

There are several functions to simplify changing population parameters, such as functions to change initial estimates and fixing parameters. As a first example, let us fix some parameters with pharmpy.modeling.fix_parameters():

PythonR

run1 = fix_parameters(model_start, ['POP_CL', 'POP_VC'])
run1.parameters

	value	lower	upper	fix
POP_CL	0.01	0.0	∞	True
POP_VC	1.00	0.0	∞	True
IIV_CL	0.10	-∞	∞	False
IIV_VC	0.10	-∞	∞	False
sigma	0.09	-∞	∞	False
IIV_CL_IIV_VC	0.01	-∞	∞	False
POP_MAT	0.10	0.0	∞	False
IIV_MAT	0.10	-∞	∞	False

Another function that may be useful would be setting the initial estimates with pharmpy.modeling.set_initial_estimates().

PythonR

run2 = set_initial_estimates(run1, {'IIV_CL': 0.05, 'IIV_VC': 0.05})
run2.parameters

	value	lower	upper	fix
POP_CL	0.01	0.0	∞	True
POP_VC	1.00	0.0	∞	True
IIV_CL	0.05	-∞	∞	False
IIV_VC	0.05	-∞	∞	False
sigma	0.09	-∞	∞	False
IIV_CL_IIV_VC	0.01	-∞	∞	False
POP_MAT	0.10	0.0	∞	False
IIV_MAT	0.10	-∞	∞	False

And then the final model code:

PythonR

print_model_code(run2)

$PROBLEM Start model
$INPUT ID TIME AMT WGT APGR DV FA1 FA2
$DATA file.csv IGNORE=@
$SUBROUTINES ADVAN2 TRANS2
$ABBR REPLACE ETA_CL=ETA(1)
$ABBR REPLACE ETA_VC=ETA(2)
$ABBR REPLACE ETA_MAT=ETA(3)
$PK
MAT = THETA(3)*EXP(ETA_MAT)
CL = THETA(1)*EXP(ETA_CL)
VC = THETA(2)*EXP(ETA_VC)
KA = 1/MAT
V = VC
$ERROR
IPRED = A(2)/VC
IF (IPRED.EQ.0) THEN
    IPREDADJ = 2.22500000000000E-16
ELSE
    IPREDADJ = IPRED
END IF
Y = IPRED + EPS(1)*IPREDADJ
$THETA  (0,0.01) FIX ; POP_CL
$THETA  (0,1.0) FIX ; POP_VC
$THETA  (0,0.1) ; POP_MAT
$OMEGA BLOCK(2)
0.05	; IIV_CL
0.01	; IIV_CL_IIV_VC
0.05	; IIV_VC
$OMEGA  0.1 ; IIV_MAT
$SIGMA  0.09 ; sigma
$TABLE ID TIME DV CIPREDI PRED CWRES FILE=mytab ONEHEADER NOAPPEND NOPRINT
$ESTIMATION METHOD=COND INTER MAXEVAL=99999

List of functions to change population parameters:

• pharmpy.modeling.add_population_parameter()

• pharmpy.modeling.fix_or_unfix_parameters()

• pharmpy.modeling.fix_parameters()

• pharmpy.modeling.fix_parameters_to()

• pharmpy.modeling.set_initial_estimates()

• pharmpy.modeling.set_lower_bounds()

• pharmpy.modeling.set_upper_bounds()

• pharmpy.modeling.unconstrain_parameters()

• pharmpy.modeling.unfix_parameters()

• pharmpy.modeling.unfix_parameters_to()

Error model

Pharmpy supports several error models. As an example, let us set the error model to a combined error model (start model had proportional error model) using pharmpy.modeling.set_combined_error_model():

PythonR

run1 = set_combined_error_model(model_start)
run1.statements.error

\begin{align*}IPRED &= \frac{A_{CENTRAL}{\left(t \right)}}{VC}\\IPREDADJ &= \begin{cases} 2.225 \cdot 10^{-16} & \text{for}\: IPRED &= 0 \\IPRED & \text{otherwise} \end{cases}\\Y &= IPRED \cdot \epsilon_{p1} + IPRED + \epsilon_{a}\\\end{align*}

List of functions to change the error model:

• pharmpy.modeling.remove_error_model()

• pharmpy.modeling.set_additive_error_model()

• pharmpy.modeling.set_combined_error_model()

• pharmpy.modeling.set_dtbs_error_model()

• pharmpy.modeling.set_iiv_on_ruv()

• pharmpy.modeling.set_power_on_ruv()

• pharmpy.modeling.set_proportional_error_model()

• pharmpy.modeling.set_time_varying_error_model()

• pharmpy.modeling.set_weighted_error_model()

• pharmpy.modeling.use_thetas_for_error_stdev()

BLQ transformations

It is also possible to perform BLQ transformations using the pharmpy.modeling.transform_blq() function. If using the M3 or M4 method the standard deviation statement is derived symbolically.

PythonR

run1 = transform_blq(model_start, method='m4', lloq=0.1)
run1.statements.error

\begin{align*}IPRED &= \frac{A_{CENTRAL}{\left(t \right)}}{VC}\\IPREDADJ &= \begin{cases} 2.225 \cdot 10^{-16} & \text{for}\: IPRED &= 0 \\IPRED & \text{otherwise} \end{cases}\\SD &= \sqrt{\sigma} \cdot \left|{IPREDADJ}\right|\\LLOQ &= 0.1\\F_{FLAG} &= \begin{cases} 0 & \text{for}\: DV \geq LLOQ \\1 & \text{otherwise} \end{cases}\\CUMD &= \phi{\left(\frac{- IPRED + LLOQ}{SD} \right)}\\CUMDZ &= \phi{\left(- \frac{IPRED}{SD} \right)}\\Y &= \begin{cases} IPRED + IPREDADJ \cdot \epsilon_{p} & \text{for}\: DV \geq LLOQ \\\frac{CUMD - CUMDZ}{1 - CUMDZ} & \text{otherwise} \end{cases}\\\end{align*}

And examine the model code:

PythonR

print_model_code(run1)

$PROBLEM Start model
$INPUT ID TIME AMT WGT APGR DV FA1 FA2
$DATA file.csv IGNORE=@
$SUBROUTINES ADVAN2 TRANS2
$ABBR REPLACE ETA_CL=ETA(1)
$ABBR REPLACE ETA_VC=ETA(2)
$ABBR REPLACE ETA_MAT=ETA(3)
$PK
MAT = THETA(3)*EXP(ETA_MAT)
CL = THETA(1)*EXP(ETA_CL)
VC = THETA(2)*EXP(ETA_VC)
KA = 1/MAT
V = VC
$ERROR
IPRED = A(2)/VC
IF (IPRED.EQ.0) THEN
    IPREDADJ = 2.22500000000000E-16
ELSE
    IPREDADJ = IPRED
END IF
SD = SQRT(SIGMA(1,1))*ABS(IPREDADJ)
LLOQ = 0.100000000000000
IF (DV.GE.LLOQ) THEN
    F_FLAG = 0
ELSE
    F_FLAG = 1
END IF
CUMD = PHI((-IPRED + LLOQ)/SD)
CUMDZ = PHI(-IPRED/SD)
IF (DV.GE.LLOQ) THEN
    Y = IPRED + EPS(1)*IPREDADJ
ELSE
    Y = (CUMD - CUMDZ)/(1 - CUMDZ)
END IF
$THETA  (0,0.01) ; POP_CL
$THETA  (0,1.0) ; POP_VC
$THETA  (0,0.1) ; POP_MAT
$OMEGA BLOCK(2)
0.1	; IIV_CL
0.01	; IIV_CL_IIV_VC
0.1	; IIV_VC
$OMEGA  0.1 ; IIV_MAT
$SIGMA  0.09 ; sigma
$TABLE ID TIME DV CIPREDI PRED CWRES FILE=mytab ONEHEADER NOAPPEND NOPRINT
$ESTIMATION METHOD=COND LAPLACE INTER MAXEVAL=99999

List of functions to perform BLQ transformations:

• pharmpy.modeling.transform_blq()

Estimation steps

Pharmpy can change the estimation steps. As an example, let us change the estimation method from FOCE to IMP and set how many iterations to output (the PRINT option in NONMEM) using the pharmpy.modeling.set_estimation_step() function:

PythonR

run1 = set_estimation_step(model_start, method='imp', keep_every_nth_iter=10)
run1.execution_steps

	method	interaction	parameter_uncertainty_method	evaluation	maximum_evaluations	laplace	isample	niter	auto	keep_every_nth_iter	tool_options
0	IMP	True	None	False	99999	False	None	None	None	10	{}

If we then examine the model code:

PythonR

print_model_code(run1)

$PROBLEM Start model
$INPUT ID TIME AMT WGT APGR DV FA1 FA2
$DATA file.csv IGNORE=@
$SUBROUTINES ADVAN2 TRANS2
$ABBR REPLACE ETA_CL=ETA(1)
$ABBR REPLACE ETA_VC=ETA(2)
$ABBR REPLACE ETA_MAT=ETA(3)
$PK
MAT = THETA(3)*EXP(ETA_MAT)
CL = THETA(1)*EXP(ETA_CL)
VC = THETA(2)*EXP(ETA_VC)
KA = 1/MAT
V = VC
$ERROR
IPRED = A(2)/VC
IF (IPRED.EQ.0) THEN
    IPREDADJ = 2.22500000000000E-16
ELSE
    IPREDADJ = IPRED
END IF
Y = IPRED + EPS(1)*IPREDADJ
$THETA  (0,0.01) ; POP_CL
$THETA  (0,1.0) ; POP_VC
$THETA  (0,0.1) ; POP_MAT
$OMEGA BLOCK(2)
0.1	; IIV_CL
0.01	; IIV_CL_IIV_VC
0.1	; IIV_VC
$OMEGA  0.1 ; IIV_MAT
$SIGMA  0.09 ; sigma
$TABLE ID TIME DV CIPREDI PRED CWRES FILE=mytab ONEHEADER NOAPPEND NOPRINT
$ESTIMATION METHOD=IMP INTER MAXEVAL=99999 PRINT=10

List of functions to change the estimation steps:

• pharmpy.modeling.add_parameter_uncertainty_step()

• pharmpy.modeling.add_estimation_step()

• pharmpy.modeling.append_estimation_step_options()

• pharmpy.modeling.remove_parameter_uncertainty_step()

• pharmpy.modeling.remove_estimation_step()

• pharmpy.modeling.set_estimation_step()

• pharmpy.modeling.set_evaluation_step()

Examining and modifying dataset

The Pharmpy dataset can be examined and modified with several help functions in Pharmpy. To read more about the dataset, see the dataset documentation.

List of dataset functions:

• pharmpy.modeling.add_time_after_dose()

• pharmpy.modeling.check_dataset()

• pharmpy.modeling.deidentify_data()

• pharmpy.modeling.drop_columns()

• pharmpy.modeling.drop_dropped_columns()

• pharmpy.modeling.expand_additional_doses()

• pharmpy.modeling.get_baselines()

• pharmpy.modeling.get_cmt()

• pharmpy.modeling.get_concentration_parameters_from_data()

• pharmpy.modeling.get_covariate_baselines()

• pharmpy.modeling.get_doseid()

• pharmpy.modeling.get_doses()

• pharmpy.modeling.get_evid()

• pharmpy.modeling.get_ids()

• pharmpy.modeling.get_mdv()

• pharmpy.modeling.get_number_of_individuals()

• pharmpy.modeling.get_number_of_observations()

• pharmpy.modeling.get_number_of_observations_per_individual()

• pharmpy.modeling.get_observations()

• pharmpy.modeling.list_time_varying_covariates()

• pharmpy.modeling.read_dataset_from_datainfo()

• pharmpy.modeling.remove_loq_data()

• pharmpy.modeling.set_covariates()

• pharmpy.modeling.set_dvid()

• pharmpy.modeling.translate_nmtran_time()

• pharmpy.modeling.undrop_columns()

Subjects

An array of all subject IDs can be retrieved.

PythonR

model = read_model(path / "pheno_real.mod")
get_ids(model)

[1,
 2,
 3,
 4,
 5,
 6,
 7,
 8,
 9,
 10,
 11,
 12,
 13,
 14,
 15,
 16,
 17,
 18,
 19,
 20,
 21,
 22,
 23,
 24,
 25,
 26,
 27,
 28,
 29,
 30,
 31,
 32,
 33,
 34,
 35,
 36,
 37,
 38,
 39,
 40,
 41,
 42,
 43,
 44,
 45,
 46,
 47,
 48,
 49,
 50,
 51,
 52,
 53,
 54,
 55,
 56,
 57,
 58,
 59]

The number of subjects in the dataset could optionally be retrieved directly.

PythonR

get_number_of_individuals(model)

59

Observations

The observations of the dataset indexed on subject ID and the independent variable can be extracted.

PythonR

get_observations(model)

ID  TIME 
1   2.0      17.3
    112.5    31.0
2   2.0       9.7
    63.5     24.6
    135.5    33.0
             ... 
58  47.5     27.9
    131.8    31.0
59  1.8      22.6
    73.8     34.3
    146.8    40.2
Name: DV, Length: 155, dtype: float64

The total number of observations can optionally be retrieved directly.

PythonR

get_number_of_observations(model)

155

Dosing

Extract dosing information

The doses of the dataset indexed on subject ID and the independent variable can be extracted.

PythonR

doses = get_doses(model)
doses

ID  TIME 
1   0.0      25.0
    12.5      3.5
    24.5      3.5
    37.0      3.5
    48.0      3.5
             ... 
59  96.0      3.0
    108.3     3.0
    120.5     3.0
    132.3     3.0
    144.8     3.0
Name: AMT, Length: 589, dtype: float64

All unique doses can be listed

PythonR

doses.unique()

array([25. ,  3.5, 15. ,  3.8, 30. ,  3.7, 18.6,  2.3, 27. ,  3.4, 24. ,
        3. , 19. ,  2.4,  3.2, 26. ,  3.3,  5. , 11. ,  2.8, 22. , 12. ,
        4. , 20. ,  2.5, 10. , 17.5,  4.5, 60. ,  7.5, 63. , 32. ,  1.9,
        1.5,  2. , 70. ,  9. , 35. , 18. , 17. ,  4.3, 34. ,  6. , 56. ,
       28. ,  7. , 14. , 16. , 40. ,  6.7,  1.7,  9.5,  4.4, 22.8])

as well as the largest and the smallest dose

PythonR

doses.min()

np.float64(1.5)

PythonR

doses.max()

np.float64(70.0)

Dose grouping

It is possible to create a DOSEID that groups each dose period starting from 1.

PythonR

ser = get_doseid(model)
ser

0       1
1       1
2       2
3       3
4       4
       ..
739    10
740    11
741    12
742    13
743    13
Name: DOSEID, Length: 744, dtype: int64

Time after dose

Add a column for time after dose (TAD)

PythonR

model = add_time_after_dose(model)
model.dataset['TAD']

0      0.0
1      2.0
2      0.0
3      0.0
4      0.0
      ... 
739    0.0
740    0.0
741    0.0
742    0.0
743    2.0
Name: TAD, Length: 744, dtype: float64

Concentration parameters

Extract pharmacokinetic concentration parameters from the dataset

PythonR

get_concentration_parameters_from_data(model)

		Cmax	Tmax	Cmin	Tmin
ID	DOSEID
1	1	17.3	2.0	NaN	NaN
	2	NaN	NaN	NaN	NaN
	3	NaN	NaN	NaN	NaN
	4	NaN	NaN	NaN	NaN
	5	NaN	NaN	NaN	NaN
...	...	...	...	...	...
59	9	NaN	NaN	NaN	NaN
	10	NaN	NaN	NaN	NaN
	11	NaN	NaN	NaN	NaN
	12	NaN	NaN	NaN	NaN
	13	40.2	2.0	NaN	NaN

589 rows × 4 columns