NONMEM PK Templates

Note

All code snippets below assume a user defined model (ADVAN9/ADVAN13 $SUBROUTINE) is applied.

Click on the green callouts to view control streams.

Absorption Models

The presence of a default dosing record in NONMEM (EVID=1 & AMT!=0) implies an instantaneous delivery of AMT of compound to the associated CMT. In some cases, the CMT associated with the dosing record happens to be the same compartment from which PK samples are taken (e.g. IV dosing with plasma PK).

However, many compounds are not delivered directly into the PK sampling compartment. Some examples include oral, subcutaneous, or inhalation administration. Conceptually, these compounds are deposited into a separate tissue space, called a DEPOT, and are then absorbed into the central (plasma) compartment with varying kinetics.

Provided below are example control streams for a variety of common absorption profiles. IV dosing will generally be described via instantaneous absorption or zero-order infusion. Other administration routes (extravascular) will likely be captured with other methods. For simplicity, all model templates in this section utilize first-order elimination from the central compartment.

Instaneous Absorption (Bolus)

While technically an IV bolus or IV push of a compound occurs over the few seconds to minutes that the clinician takes to depress the syringe or otherwise administer the drug, modeling the drug absorption as instantaneous can be useful when the timescale of delivery is short enough to not influence the PK sampling results, or if there is insufficient data to support the added complexity of modeling an infusion. If the IV delivery takes place over a relatively longer duration which may influence the PK results, an infusion model (described under ‘Zero-Order Absorption’) is likely more appropriate.

Instantaneous Bolus into Central Compartment

The RATE column must be missing from the dataset or set to 0 on dosing records. The dosing records and PK records will be associated with the same CMT (in this case, 1)

Note that the \(\$\)DES block contains no information about the input. Dosing events are detailed in the dataset.

$PROBLEM Bolus Input into Central Compartment with First-Order Elimination
;Bolus into Sampling CMT
;First-Order Elimination (k10)
;CENTRAL=CMT1 (No Depot)

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(CENTRAL)
 
$PK
;Typical Values
TVVC = THETA(1)
TVCL = THETA(2)

;Parameters
VC = TVVC*EXP(ETA(1))
CL = TVCL*EXP(ETA(2))

;Derived Parameters
K10 = CL/VC

$DES
DADT(1) = -K10*A(1)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Zero-Order Absorption (Infusions)

Zero-order absorption models can be conceptualized as infusions. These infusions can be delivered directly to the central compartment (IV administration), or infused into a depot compartment that is subsequently absorbed into the central compartment.

For infusions, \(infusion~rate = \frac{amount}{infusion~duration}\). To fully describe the dosing event, either the rate (R) or duration (D) must be defined in the control stream as R(N) or D(N), respectively, where (N) is replaced with the CMT receiving the infusion.

Depending on the information available to the PK analyst, the rate or duration can either be estimated as a parameter (THETA) or provided as a column in the NONMEM dataset. For example, when compounds are delivered via IV, the PK analyst frequently has access to the exact durations of the IV infusions, which can directly be assigned to D (e.g. D1 to define the duration into CMT1), reducing the number of model-estimated parameters.

The RATE (rate of infusion) column enables the analyst to decide if either the rate or duration will be provided/estimated. On dosing records, the allowed values of RATE are 0, -1, and -2.

RATE	Description
0	Bolus dose, not an infusion
-1	Infusion dose, infusion rate (R) estimated/provided in control stream
-2	Infusion dose, infusion duration (D) estimated/provided in control stream
.	Rate must be missing on other/observation records

Regardless of the analyst choice of RATE, one can convert between infusion rate and infusion duration from model outputs if the AMT is known. Examples of zero-order absorption control streams are provided below.

Zero-Order Infusion Into Central Compartment (IV)

When the infusion is modeled as directly entering the central compartment, a “depot” does not have to be explicitly modeled as a separate compartment. Note that the CMT of the PK observations (in this case, 1) will be the same as the CMT associated with the dosing records.

For these examples, note that the \(\$\)DES block has no input rate. Dosing inputs are handled in the dataset and via declaration of D or R.

Infusion Rates Estimated

RATE must be set to -1 on dosing records. R1 can be declared with or without inter-individual variability. In this example, there is IIV present on R1.

$PROBLEM Zero-Order Infusion (R1 Estimated) into Central Compartment
;Zero-Order Absorption Rate (R1)
;First-Order Elimination (k10)
;CENTRAL=CMT1 (No Depot)

$INPUT IGNORE=C ID TIME DV AMT EVID CMT RATE

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(CENTRAL)
 
$PK
;Typical Values
TVR1 = THETA(1)
TVVC = THETA(2)
TVCL = THETA(3)

;Parameters
R1 = TVR1*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
CL = TVCL*EXP(ETA(3))

;Derived Parameters
K10 = CL/VC

$DES
DADT(1) = -K10*A(1)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Infusion Rates Provided

RATE must be set to -1 on dosing records.

$PROBLEM Zero-Order Infusion (R1 Provided) into Central Compartment
;Zero-Order Absorption Rate (R1)
;First-Order Elimination (k10)
;CENTRAL=CMT1 (No Depot)
;Infusion rates (units of mass per time) provided in INFRATE

$INPUT IGNORE=C ID TIME DV AMT EVID CMT RATE INFRATE

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9
$MODEL
COMP=(CENTRAL)
 
$PK
;Typical Values
TVVC = THETA(1)
TVCL = THETA(2)

;Parameters
R1 = INFRATE
VC = TVVC*EXP(ETA(1))
CL = TVCL*EXP(ETA(2))

;Derived Parameters
K10 = CL/VC

$DES
DADT(1) = -K10*A(1)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Infusion Durations Estimated

RATE must be set to -2 on dosing records. D1 can be declared with or without inter-individual variability. In this example, there is IIV present on D1.

$PROBLEM Zero-Order Infusion (D1 Estimated) into Central Compartment
;Duration of Zero-Order Infusion (D1)
;First-Order Elimination (k10)
;CENTRAL=CMT1 (No Depot)

$INPUT IGNORE=C ID TIME DV AMT EVID CMT RATE

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(CENTRAL)
 
$PK
;Typical Values
TVD1 = THETA(1)
TVVC = THETA(2)
TVCL = THETA(3)

;Parameters
D1 = TVD1*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
CL = TVCL*EXP(ETA(3))

;Derived Parameters
K10 = CL/VC

$DES
DADT(1) = -K10*A(1)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Infusion Durations Provided

RATE must be set to -2 on dosing records.

$PROBLEM Zero-Order Infusion (D1 Provided) into Central Compartment
;Duration of Zero-Order Infusion (D1)
;First-Order Elimination (k10)
;CENTRAL=CMT1 (No Depot)
;Infusion durations (units of time) provided in DUR

$INPUT IGNORE=C ID TIME DV AMT EVID CMT RATE DUR

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(CENTRAL)
 
$PK
;Typical Values
TVVC = THETA(1)
TVCL = THETA(2)

;Parameters
D1 = DUR
VC = TVVC*EXP(ETA(1))
CL = TVCL*EXP(ETA(2))

;Derived Parameters
K10 = CL/VC

$DES
DADT(1) = -K10*A(1)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Zero-Order Infusion into Depot

This type of model is useful for plasma PK profiles that exhibit sigmoidal character during absorption, or when a first-order + lag time model presents with numerical instability. In the examples below, the depot compartment receives the zero-order infusion. Drug transfer from the depot to the central compartment occurs through first-order absorption.

Conceptually, this may be useful for subcutaneous or subdermal compounds that, intially, slowly leach into plasma but after a delay are absorbed more quickly, perhaps after some dissolution or dispersion occurs.

It is unlikely that any infusion rates/durations available from clinical data will directly translate to parameters in the model structure, as the conceptual depot compartment will not necessarily reflect the anatomical depot location. As such, it is preferable to estimate R1 or D1 as a fixed effect parameter within the control stream.

Infusion Rates Estimated

RATE must be set to -1 on dosing records. R1 can be declared with or without inter-individual variability. In this example, there is IIV present on R1.

$PROBLEM Zero-Order Infusion (R1 Estimated) into Depot Compartment
;Zero-Order Infusion Rate into Depot (R1)
;First-Order Absorption from Depot into Central Compartment (k12)
;First-Order Elimination from Central Compartment (k20)
;DEPOT=CMT1, CENTRAL=CMT2
;Note that dosing records will have CMT=1 and the PK records will have CMT=2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT RATE

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Typical Values
TVR1 = THETA(1)
TVKA = THETA(2)
TVVC = THETA(3)
TVCL = THETA(4)

;Parameters
R1 = TVR1*EXP(ETA(1))
KA = TVKA*EXP(ETA(2))
VC = TVVC*EXP(ETA(3))
CL = TVCL*EXP(ETA(4))

;Derived Parameters
K12 = KA
K20 = CL/VC

$DES
DADT(1) = -K12*A(1)
DADT(2) = K12*A(1)-K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Infusion Durations Estimated

RATE must be set to -2 on dosing records. D1 can be declared with or without inter-individual variability. In this example, there is IIV present on D1.

$PROBLEM Zero-Order Infusion (D1 Estimated) into Depot Compartment
;Zero-Order Infusion Duration into Depot (D1)
;First-Order Absorption from Depot into Central Compartment (k12)
;First-Order Elimination from Central Compartment (k20)
;DEPOT=CMT1, CENTRAL=CMT2
;Note that dosing records will have CMT=1 and the PK records will have CMT=2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT RATE

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Typical Values
TVD1 = THETA(1)
TVKA = THETA(2)
TVVC = THETA(3)
TVCL = THETA(4)

;Parameters
D1 = TVD1*EXP(ETA(1))
KA = TVKA*EXP(ETA(2))
VC = TVVC*EXP(ETA(3))
CL = TVCL*EXP(ETA(4))

;Derived Parameters
K12 = KA
K20 = CL/VC

$DES
DADT(1) = -K12*A(1)
DADT(2) = K12*A(1)-K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

First-Order Absorption

First-order absorption models are often the primary starting point for extravascularly administered compounds.

Per Fick’s law of diffusion, mass moves from regions of high concentration (i.e. DEPOT) to low concentration (i.e. CENTRAL) through permeable membranes at a rate proportional to the concentration difference \(\Delta{}C = (C_{central} - C_{depot})\). Assuming the concentration in the CENTRAL compartment \(C_{central}\) is negligible compared to the concentration in the DEPOT compartment \(C_{depot}\) such that \((C_{central} - C_{depot}) ~\sim~-C_{depot}\), and that no back-transfer from CENTRAL to the DEPOT can occur, the concentration of drug in the depot compartment over time simplifies to \(\frac{dC}{dt} = -k(C_{depot})\), where k is a function of the permeability between the compartments.

When modeling, the volume of the depot compartment is ignored, and the expression is finally simplified to \(\frac{d(A1)}{dt} = -k(A1)\), where \(A1\) refers to the mass in the depot compartment. Through conservation of mass, the corresponding expression for the input into the central compartment over time is \(\frac{d(A2)}{dt} = k(A1)\)

First-Order Absorption

The dosing records will be associated with the DEPOT compartment, in this case, CMT1. The PK records will be associated with the CENTRAL compartment, in this case, CMT2.

At the TIME of dosing records, an instantaneous amount AMT will be added to the DEPOT compartment.

;First-Order Absorption (ka) from Depot into Central
;First-Order Elimination (ke) from Central
;DEPOT=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT RATE

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Typical Values
TVKA = THETA(1)
TVVC = THETA(2)
TVCL = THETA(3)

;Parameters
KA = TVKA*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
CL = TVCL*EXP(ETA(3))

;Derived Parameters
K12 = KA
K20 = CL/VC

$DES
DADT(1) = -K12*A(1)
DADT(2) =  K12*A(1) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Lag Times

Often, there is a notable delay between the TIME associated with dosing and the time at which plasma concentrations begin to rise. An example would be an orally delivered drug primarily transported through the lower intestine, which after dosing (swallowing) would require passage through the GI tract before being absorbed in the body. A common approach to model this delay is through a LAG time.

In practice, associating a lag time to a compartment simply adds a delay (ALAG) between the TIME present on the dosing record and the time at which the AMT appears within the corresponding compartment.

Lag times can be added to compartments by declaring a value for ALAG(N) within the control stream, where (N) is replaced with the CMT number associated with the delay. Lag times are virtually always estimated as a fixed effect parameter (THETA).

Lag time models, especially if including IIV on ALAG, frequently result in numerical difficulties given the abrupt on/off nature of the central compartment input. Alternatives to the lag time model include transit compartment absorption and zero-order infusion into a depot. However, there is tremendous value in the basic lag-time approach.

1. There is only one parameter to estimate.
2. ALAG is directly interpretable.
3. It is popular and easily communicated to clinical colleagues/reviewers.

The value of simplicity should be a factor when deciding between a LAG and a more complicated model structure (in addition to statistical fit).

First-Order Absorption With Lag Time

The dosing records will be associated with the DEPOT compartment, in this case, CMT1. The PK records will be associated with the CENTRAL compartment, in this case, CMT2.

At the TIME of dosing records + ALAG1, an amount AMT will be added to the DEPOT compartment.

IIV can be added on ALAG. In this case, there is no IIV on ALAG.

;First-Order Absorption (ka) from Depot into Central with Lag time
;First-Order Elimination (ke) from Central
;DEPOT=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Typical Values
TVKA = THETA(1)
TVALAG1 = THETA(2)
TVVC = THETA(3)
TVCL = THETA(4)

;Parameters
KA = TVKA*EXP(ETA(1))
ALAG1 = TVALAG1
VC = TVVC*EXP(ETA(2))
CL = TVCL*EXP(ETA(3))

;Derived Parameters
K12 = KA
K20 = CL/VC

$DES
DADT(1) = -K12*A(1)
DADT(2) =  K12*A(1) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Dual Absorption (Zero-Order/First-Order/LAG)

While simple zero-order infusions and first-order absorption models frequently provide adequate fits to pharmacokinetic data, sometimes the data indicate that a mixture of two absorption models, with or without a delay, are necessary. One obvious sign suggesting a multiple absorption model may be necessary is the presence of multiple ‘peaks’ in the concentration-time profile. An imbalanced distribution of residuals specific to parts of the absorption phase may also indicate a need to adopt a mixed absorption model.

For oral compounds, mixed absorption can be conceptualized as reflecting multiple modes of transport, like passive diffusion and transporter-mediated influx. It may also reflect changes in the ability of the drug to enter the body along different sections of the GI tract, due to pH, the gradient of intestinal transporter density, or microbial enzymatic processes, etc. For subcutaneous administration, dual absorption may reflect a ‘partitioning’ of the drug within anatomical sections of the subcutaneous space which differentially absorb into plasma, or potential chemical/physical changes in certain fractions of the deposited dose, like aggregation of a protein biotherapeutic or crystallization of a small molecule after diffusion of the vehicle.

Practically, encoding a dual absorption model will require a NONMEM-ready dataset with two dosing records at each actual clinical dosing event, one corresponding to each CMT where drug is to be administered. Here, examples of dosing records for dual absorption are provided for models with one or two depot compartments, where the subject receives a single 100 mg dose at TIME=0.

Example Dosing Record, Modeling Two Depots

Here, CMT1=DEPOT1, CMT2=DEPOT2, and CMT3=PLASMA PK. Part of the 100 mg dose is administered to CMT1 and the remainder is administered to CMT2. The dose can be absorbed from CMT1 and CMT2 into the plasma compartment, CMT3, with differing kinetics or lag times.

C	ID	TIME	DV	AMT	EVID	CMT
.	1	0	.	100	1	1
.	1	0	.	100	1	2
.	1	0.25	1.45	.	0	3
.	1	0.5	6.88	.	0	3
.	1	1	20.57	.	0	3

Example Dosing Record, Modeling One Depot and One Direct Infusion to Plasma

Here, part of the dose will be administered to a DEPOT (CMT1), and another part of the dose will be administered directly to the CENTRAL/PLASMA PK CMT (CMT2). For more on administration directly to the plasma compartment, see the earlier sections on Bolus and Infusion absorption models.

^a0 for IV Bolus, -1 when specifying infusion rate, -2 when specifying infusion duration

C	ID	TIME	DV	AMT	EVID	CMT	RATE
.	1	0	.	100	1	1	0
.	1	0	.	100	1	2	-2a
.	1	0.25	1.45	.	0	2	.
.	1	0.5	6.88	.	0	2	.
.	1	1	20.57	.	0	2	.

At first glance, it seems the subject receives a full dose into each specified CMT, thus receiving twice as much drug as appropriate. However, the correct partitioning of the dose is taken care of in the control stream. A fraction of the dose, F(N), is defined for each CMT receiving drug, where (N) is replaced with the corresponding CMT. For example, the fraction of AMT allocated to CMT1 is designated F1, and the fraction of AMT allocated to CMT2 is designated F2.

F2 can be defined in terms of F1, \(F2 = 1 - F1\), since \(F1 + F2 = 1\). In this way, although the dosing record appears twice, the sum of the deposited AMTs into each DEPOT correctly reflect the clinical AMT administered at TIME=0.

F1 can be estimated as a fixed effect parameter THETA. The estimation boundaries of F1 are (0, 1).

As values of F(N) are bounded between (0, 1), inclusion of inter-individual variability on F1 necessitates a logit-transformation to ensure individual estimates of F(N) remain within (0, 1).

Two First-Order Absorptions Beginning Simultaneously at Dosing

For each clinical dose, there will be two dosing records, once for each DEPOT compartment, in this case, CMT1 and CMT2. The PK records will be associated with the CENTRAL compartment, in this case, CMT3.

In this model, fractions of the dose F1 and F2 are absorbed via first-order absorption rate constants ka1 and ka2, respectively, both beginning at the time of dosing.

$PROBLEM Dual First-Order Absorption Both Starting At Dosing
;First-Order Absorption (ka1) from Depot CMT1 starting at dosing
;First-Order Absorption (ka2) from Depot CMT2 starting at dosing
;First-Order Elimination (ke) from Central
;DEPOT1=CMT1, DEPOT2=CMT2, CENTRAL=CMT3

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9


$MODEL
COMP=(DEPOT1)
COMP=(DEPOT2)
COMP=(CENTRAL)
 
$PK
;Typical Values
TVKA1 = THETA(1)
TVKA2 = THETA(2)
TVVC = THETA(3)
TVCL = THETA(4)
TVF1 = THETA(5)

;Parameters
KA1 = TVKA1*EXP(ETA(1))
KA2 = TVKA2*EXP(ETA(2))
VC = TVVC*EXP(ETA(3))
CL = TVCL*EXP(ETA(4))

F1 = TVF1
F2 = 1 - F1

;Derived Parameters
K13 = KA1
K23 = KA2
K30 = CL/VC

$DES
DADT(1) = -K13*A(1)
DADT(2) = -K23*A(2)
DADT(3) =  K13*A(1) + K23*A(2) - K30*A(3)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Dual First-Order and Zero-Order Absorptions Beginning Simultaneously at Dosing

For each clinical dose, there will be two dosing records, once for each compartment receiving drug, in this case, CMT1 and CMT2. The PK records will be associated with the CENTRAL compartment, in this case, CMT2.

In this model, a fraction of the dose F1 is absorbed via first-order absorption with rate constant ka1. The remaining fraction of the dose, F2, undergoes zero-order infusion directly into the plasma compartment CMT2 with estimated duration D2.

;First-Order Absorption (ka1) from Depot CMT1 starting at dosing
;Zero-order infusion with duration D2 into CMT2 (DEPOT2 and PK) starting after dosing
;First-Order Elimination (ke) from Central
;DEPOT1=CMT1, DEPOT2=CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT RATE

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Typical Values
TVKA = THETA(1)
TVD2 = THETA(2)
TVVC = THETA(3)
TVCL = THETA(4)
TVF1 = THETA(5)

;Parameters
KA = TVKA*EXP(ETA(1))
D2 = TVD2
VC = TVVC*EXP(ETA(2))
CL = TVCL*EXP(ETA(3))

F1 = TVF1
F2 = 1 - F1

;Derived Parameters
K12 = KA
K20 = CL/VC

$DES
DADT(1) = -K12*A(1)
DADT(2) =  K12*A(1) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

The dual absorptions do not have to begin simultaneously. Specific fractions can undergo delayed absorption by adding lag times on the target CMT.

Two First-Order Absorptions, One Beginning After a Delay

For each clinical dose, there will be two dosing records, once for each DEPOT compartment, in this case, CMT1 and CMT2. The PK records will be associated with the CENTRAL compartment, in this case, CMT3.

In this model, a fraction of the dose F1 is absorbed via first-order absorption beginning at dosing. Another fraction of the dose F2 is absorbed via first-order absorption after a lag time ALAG2. At the TIME of dosing records, an instantaneous amount AMT*F1 will be added to DEPOT compartment CMT1. At the TIME of dosing records + ALAG2, an amount AMT*F2 will be added to DEPOT compartment CMT2.

;First-Order Absorption (ka1) from Depot CMT1 starting at dosing
;First-Order Absorption (ka2) from Depot CMT2 starting after ALAG2
;First-Order Elimination (ke) from Central
;DEPOT1=CMT1, DEPOT2=CMT2, CENTRAL=CMT3

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT1)
COMP=(DEPOT2)
COMP=(CENTRAL)
 
$PK
;Typical Values
TVKA1 = THETA(1)
TVKA2 = THETA(2)
TVALAG2 = THETA(3)
TVVC = THETA(4)
TVCL = THETA(5)
TVF1 = THETA(6)

;Parameters
KA1 = TVKA1*EXP(ETA(1))
KA2 = TVKA2*EXP(ETA(2))
ALAG2 = TVALAG2
VC = TVVC*EXP(ETA(3))
CL = TVCL*EXP(ETA(4))

F1 = TVF1
F2 = 1 - F1

;Derived Parameters
K13 = KA1
K23 = KA2
K30 = CL/VC

$DES
DADT(1) = -K13*A(1)
DADT(2) = -K23*A(2)
DADT(3) =  K13*A(1) + K23*A(2) - K30*A(3)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Dual First-Order and Zero-Order Absorption, with the Zero-Order Process Delayed

For each clinical dose, there will be two dosing records, once for each DEPOT compartment, in this case, CMT1 and CMT2. The RATE column for dosing occasions associated with CMT1 must be equal to 0. Since we are estimating the duration of infusion D2, the RATE column for dosing occasions associated with CMT2 must be set to -2. The PK records will be associated with the CENTRAL compartment, in this case, CMT2 (same as the 2nd depot).

In this model, a fraction of the dose F1 is absorbed via first-order absorption after dosing. Another fraction of the dose F2 is delivered as a zero-order infusion after a lag time ALAG2. At the TIME of dosing records, an instantaneous amount AMT*F1 will be added to DEPOT compartment CMT1. At the TIME of dosing records + ALAG2, a zero-order infusion of AMT*F2 with duration D2 begins. Here, no IIV is estimated on D2.

;First-Order Absorption (ka1) from Depot CMT1 starting at dosing
;Zero-order infusion with duration D2 into CMT2 (DEPOT2 and PK) starting after ALAG2
;First-Order Elimination (ke) from Central
;DEPOT1=CMT1, DEPOT2=CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT RATE

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Typical Values
TVKA = THETA(1)
TVD2 = THETA(2)
TVALAG2 = THETA(3)
TVVC = THETA(4)
TVCL = THETA(5)
TVF1 = THETA(6)

;Parameters
KA = TVKA*EXP(ETA(1))
D2 = TVD2
ALAG2 = TVALAG2
VC = TVVC*EXP(ETA(2))
CL = TVCL*EXP(ETA(3))

F1 = TVF1
F2 = 1 - F1

;Derived Parameters
K12 = KA
K20 = CL/VC

$DES
DADT(1) = -K12*A(1)
DADT(2) =  K12*A(1) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Dual First-Order and Zero-Order Absorption, with the First-Order Process Delayed

This is identical to the previous control stream except that the lag time appears on the first-order depot compartment CMT1.

;First-Order Absorption (ka) from Depot CMT1 starting after ALAG1
;Zero-order infusion with duration D2 into CMT2 (DEPOT2 and PK)
;First-Order Elimination (ke) from Central
;DEPOT1=CMT1, DEPOT2=CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT RATE

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Typical Values
TVKA = THETA(1)
TVD2 = THETA(2)
TVALAG1 = THETA(3)
TVVC = THETA(4)
TVCL = THETA(5)
TVF1 = THETA(6)

;Parameters
KA = TVKA*EXP(ETA(1))
D2 = TVD2
ALAG1 = TVALAG1
VC = TVVC*EXP(ETA(2))
CL = TVCL*EXP(ETA(3))

F1 = TVF1
F2 = 1 - F1

;Derived Parameters
K12 = KA
K20 = CL/VC

$DES
DADT(1) = -K12*A(1)
DADT(2) =  K12*A(1) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Including IIV on F1 (Logit-transformation)

For any of the models above, IIV can be included on F1 through a logit-transformation. The code snippet below can be included in the \(\$\)PK block, replacing the previous assignments of F1 and F2.

Here, TVF1, the typical value of the fraction of dose allocated to DEPOT1, is still constrained between 0 and 1. LGTF1 can vary between -Inf and +Inf. F1 can vary between individuals and is constrained to be between 0 and 1. Note that the THETA associated with TVF1 is on the identity scale, but the IIV is normally distributed on the logit-scale.

$PK
TVF1 = THETA(N) ;Substitute (N) as appropriate

;Logit Calculations
PHI = TVF1 / (1 - TVF1)
LGTF1 = LOG(PHI) + ETA(N) ;Substitute (N) as appropriate

F1 = EXP(LGTF1) / (1 + EXP(LGTF1))
F2 = 1 - F1

Sequential Zero-Order then First-Order Absorption (Non-Overlapping)

The sources of absorption can be set to start one after the other. For the control streams below, the moment a zero order infusion ends, a first order absorption begins. This is achieved by setting the lag time of the first-order absorption compartment equal to the duration of the zero-order infusion D2.

Sequential Switch from Zero-Order to a First-Order Absorption

;First-Order Absorption (ka) from Depot CMT1 starting after ALAG1
;Zero-order infusion with duration D2, equal to ALAG1, into CMT2
;First-Order Elimination (ke) from Central
;DEPOT1=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT RATE

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Typical Values
TVKA = THETA(1)
TVD2 = THETA(2)
TVVC = THETA(3)
TVCL = THETA(4)
TVF1 = THETA(5)

;Parameters
KA = TVKA*EXP(ETA(1))
D2 = TVD2
ALAG1 = D2 ;The delay before first order is equal to the duration of the zero-order process 
VC = TVVC*EXP(ETA(2))
CL = TVCL*EXP(ETA(3))

F1 = TVF1 ;Fraction absorbed via first-order
F2 = 1 - F1 ;Fraction absorbed via zero-order

;Derived Parameters
K12 = KA
K20 = CL/VC

$DES
DADT(1) = -K12*A(1)
DADT(2) =  K12*A(1) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Linked Sequential Switch from Zero-Order to a First-Order Absorption

Here, the absorptions are ‘Linked’ by constraining the first-order absorption rate at the time of the zero-order/first-order switch to be equal to the rate throughout the zero-order absorption. This trick can be useful for oral compounds which exhibit solubility limited absorption (zero-order), and once fully solubilized are taken up through first-absorption from the same depot compartment (gut). See Holford et al. (1992) J. Pharmacokinetics and Biopharmaceutics for more.

A benefit of this model compared to the independent sequential switch is that the first order absorption rate constant, ka, is a function of F1 and D2, and does not have to be independently estimated.

Here, the zero-order process starts at TIME 0. If a delay on the zero-order infusion (ALAG2) is required, note that the lag time on the first-order absorption ALAG1 would then be set as ALAG1 = D2 + ALAG2.

;First-Order Absorption (ka) from Depot CMT1 starting after ALAG1
;ka is defined as a function of F1 and D2
;Zero-order infusion with duration D2, equal to ALAG1, into CMT2
;First-Order Elimination (ke) from Central
;DEPOT1=CMT1, DEPOT2=CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT RATE

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Typical Values
TVD2 = THETA(1)
TVVC = THETA(2)
TVCL = THETA(3)
TVF1 = THETA(4)

;Parameters
D2 = TVD2 
ALAG1 = D2 ;The delay before first order is equal to the duration of the zero-order process 
VC = TVVC*EXP(ETA(1))
CL = TVCL*EXP(ETA(2))

F1 = TVF1 ;Fraction absorbed by first order
F2 = 1 - F1 ;Fraction absorbed by zero-order

KA = F2/(F1 * D2)  

;Derived Parameters
K12 = KA
K20 = CL/VC

$DES
DADT(1) = -K12*A(1)
DADT(2) =  K12*A(1) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Saturable (Michaelis-Menten) Absorption

Saturable Michaelis-Menten Absorption

The dosing records will be associated with the DEPOT compartment, in this case, CMT1. The PK records will be associated with the CENTRAL compartment, in this case, CMT2.

;Michaelis-Menten Absorption VMAX/KM from Depot into Central
;First-Order Elimination (ke) from Central
;DEPOT=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Typical Values
TVVMAX = THETA(1)
TVKM = THETA(2)
TVVC = THETA(2)
TVCL = THETA(3)

;Parameters
VMAX = TVVMAX*EXP(ETA(1))
KM = TVKM*EXP(ETA(2))
VC = TVVC*EXP(ETA(3))
CL = TVCL*EXP(ETA(4))

;Derived Parameters
K20 = CL/VC

$DES
DADT(1) = -VMAX*A(1)/(KM + A(1))
DADT(2) =  VMAX*A(1)/(KM + A(1)) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Transit Compartment Models

A popular alternative to first order models with lag times are transit compartment models. As one may suspect, transit compartment models reflect a series of linked compartments that the drug passes through prior to being delivered to the CENTRAL compartment. The absorption kinetics are then defined by the number of transit compartments \(N_{cpt}\), the mean transit time through the transit system \(MTT\), and the transfer rate between compartment \(k_{tr}\). Fortunately, the transit rate can be often be defined in terms of \(N_{cpt}\) and \(MTT\).

There are inconsistencies in the pharmacometric literature about the definition of \(N_{cpt}\). For the purposes of this reference, \(N_{cpt}\) will correspond to the number of transit compartments that receive both first-order input and undergo first-order elimination, i.e. \(\frac{dA(n)}{dt}=k_{tr}*A(n-1) - k_{tr}*A(n)\).

The first “compartment” in the transit chain, which receives a bolus dose and only undergoes first-order elimination i.e. \(\frac{dA_{0}}{dt}=-k_{tr}*A(0)\) DOES NOT COUNT TOWARDS \(N_{cpt}\).

Transit models have certain benefits over lag time models. First, they can capture a slow trickle of absorption that would occur during a potential lag time component, and second, they avoid discontinuties in input rates that can lead to numerical difficulties, especially when considering IIV. The downside is they require a bit more optimization and analytical solutions require certain assumptions to be met if multiple doses are administered.

Erlang Transit Distribution (Analytical)

The Erlang transit model is the analytical solution to the input rate into the central compartment of a series of transit compartments where the number of compartments is set to a predefined integer value. Multiple models with differing integer numbers of compartments can be run and compared for statistical fit.

As the analytical function describes the input into the central compartment, the DEPOT compartment is completely ignored. The control stream presented below encodes a model in which the linked transit compartments deposit directly into the central compartment (as opposed to depositing into a DEPOT which would then undergo first order absorption with rate constant \(k_a\) where \(k_a \ne k_{tr}\)). This control stream will work for multiple doses, with the assumption that virtually all of the drug has been absorbed into the central compartment prior to the next dose.

For the Erlang model, \(k_{tr} = \frac{N_{cpt}+1}{MTT}\), where \(N_{cpt}\) does not include the initial compartment receiving the bolus dose.

\(f(t) = Dose*\frac{k_{tr}^{N_{cpt} + 1}}{\Gamma{}(N_{cpt} + 1)}*t^{N_{cpt}}*e^{-k_{tr}t}\), where \(f(t)\) is the instantaneous rate of amount exiting the final transit compartment \(\frac{-dA(N_{cpt})}{dt}\).

The value of the gamma function

\(\Gamma{}(N) = (N - 1)!\)

We can make use of the GAMLN function in NONMEM, which returns the natural logarithm of the gamma function of the input.

The dosing records will be associated with the DEPOT compartment, in this case, CMT1. The PK records will be associated with the CENTRAL compartment, in this case, CMT2.

;Erlang Distibution for 4 Transit Compartments (Not including the initial compartment receiving the bolus dose)
;First-Order Elimination (ke) from Central
;DEPOT=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Update Doses and Dose Times
IF(AMT.GT.0)DOSE = AMT
IF(AMT.GT.0)TDOSE = TIME

;Erlang Inputs
NCPT = 4 ;Set to number of compartments
TVMTT = THETA(1) ;Mean Transit Time
F1 = 0; The value in the depot compartment does not matter

;Typical Values
TVVC = THETA(2)
TVCL = THETA(3)

;Parameters
MTT = TVMTT*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
CL = TVCL*EXP(ETA(3))

;Derived Parameters
KTR = (NCPT + 1)/MTT ;Transit rate through the linked compartments
K20 = CL/VC
GAMNCPT = EXP(GAMLN(NCPT + 1)) ;Gamma function of the number of transit compartments + 1

$DES
TCUR = T - TDOSE ;Time after dosing

DADT(1) = 0
DADT(2) = DOSE*(KTR**(NCPT + 1)/GAMNCPT)*(TCUR**(NCPT))*EXP(-KTR*TCUR) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Erlang Transit Distribution (Analytical - Into a Separate Absorption Compartment)

This transit model is the analytical solution to the input rate into a distinct absorption compartment from a series of transit compartments where the number of compartments is set to a predefined integer value. The distinct absorption compartment then undergoes first order absorption into the central compartment with an absorption rate \(k_a \ne k_{tr}\) which is separately estimated from the transit rate constant \(k_{tr}\). Multiple models with differing integer numbers of compartments can be run and compared for statistical fit. This control stream will work for multiple doses, with the assumption that virtually all of the drug has been absorbed into the central compartment prior to the next dose.

For the Erlang model, \(k_{tr} = \frac{N_{cpt}+1}{MTT}\), where \(N_{cpt}\) does not include the initial compartment receiving the bolus dose.

\(f(t) = Dose*\frac{k_{tr}^{N_{cpt} + 1}}{\Gamma{}(N_{cpt} + 1)}*t^{N_{cpt}}*e^{-k_{tr}t}\), where \(f(t)\) is the instantaneous rate of amount exiting the final transit compartment \(\frac{-dA(N_{cpt})}{dt}\).

The value of the gamma function

\(\Gamma{}(N) = (N - 1)!\)

We can make use of the GAMLN function in NONMEM, which returns the natural logarithm of the gamma function of the input.

The dosing records will be associated with the DEPOT compartment, in this case, CMT1. The PK records will be associated with the CENTRAL compartment, in this case, CMT2.

;Erlang Distibution for 4 Transit Compartments (not including the first compartment which receives the bolus dose or the DEPOT which enters the CENTRAL)
;First-Order Elimination (ke) from Central
;DEPOT=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Update Doses and Dose Times
IF(AMT.GT.0)DOSE = AMT
IF(AMT.GT.0)TDOSE = TIME

;Erlang Inputs
NCPT = 4 ;Set to number of compartments
TVMTT = THETA(1) ;Mean Transit Time
F1 = 0; We do not use the AMT in the dosing record as the input is analytically derived
KA = THETA(2)

;Typical Values
TVVC = THETA(3)
TVCL = THETA(4)

;Parameters
MTT = TVMTT*EXP(ETA(1))
KA = TVKA*EXP(ETA(2))
VC = TVVC*EXP(ETA(3))
CL = TVCL*EXP(ETA(4))

;Derived Parameters
KTR = (NCPT+1)/MTT ;Transit rate through the linked compartments
K12 = KA ;absorption rate constant from DEPOT to CENTRAL
K20 = CL/VC
GAMNCMT = EXP(GAMLN(NCPT + 1)) ;Gamma function of the number of transit compartments + 1

$DES
TCUR = T - TDOSE ;Time after dosing
DADT(1) = DOSE*(KTR**(NCPT + 1)/GAMNCPT)*(TCUR**(NCPT))*EXP(-KTR*TCUR) - K12*A(1)
DADT(2) = K12*A(1) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Transit Compartment Modeling (Compartments Explicitly Defined)

This model explicitly defines the number of transit compartments, which may be necessary in the case of multiple dosing (where the assumption of full absorption prior to dosing is not met). The control stream below is an example for which 4 transit compartments are defined.

The dosing records will be associated with the compartment where the dose is initialized, in this case, CMT1. This compartment does not count towards \(N_{cpt}\). The PK records will be associated with the CENTRAL compartment, in this case, CMT2.

CMTs 3, 4, 5, and 6, represent the transit compartment chain (4 compartments total). CMT1 will feed into CMT3, which feeds to 4, which feeds to 5, which feeds to 6, which then feeds into the central compartment CMT2.

;4 Transit Compartments into Central
;First-Order Elimination (ke) from Central
;Transit Compartment 0=CMT1, CENTRAL=CMT2
;Transit Compartments 1, 2, 3, 4 = CMT3/4/5/6 

$INPUT IGNORE=C ID TIME DV AMT EVID CMT RATE

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DOSING)
COMP=(CENTRAL)
COMP=(TRANS1)
COMP=(TRANS2)
COMP=(TRANS3)
COMP=(TRANS4)
 
$PK
NCPT=4 ;4 Transit compartments
TVMTT = THETA(1) ;Mean Transit Time

;Typical Values
TVVC = THETA(2)
TVCL = THETA(3)

;Parameters
MTT = TVMTT*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
CL = TVCL*EXP(ETA(3))

;Derived Parameters
KTR = (NCPT+1)/MTT ;Transit rate through the linked compartments
K20 = CL/VC

$DES
DADT(1) = -KTR*A(1) ;Dosing compartment
DADT(3) = KTR*(A(1) - A(3)) ;Transit compartment 1
DADT(4) = KTR*(A(3) - A(4)) ;Transit compartment 2
DADT(5) = KTR*(A(4) - A(5)) ;Transit compartment 3
DADT(6) = KTR*(A(5) - A(6)) ;Transit compartment 4
DADT(2) = KTR*A(6) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Transit Compartment Modeling (Compartments Explicitly Defined - Into a Distinct Absorption Compartment)

This model explicitly defines the number of transit compartments, which may be necessary in the case of multiple dosing (where the assumption of full absorption prior to dosing is not met). The control stream below is an example for which 4 transit compartments are defined. The final transit compartment deposits into a distinct absorption compartment which then undergoes first order absorption into the central compartment with an absorption rate \(k_a \ne k_{tr}\) which is separately estimated from the transit rate constant \(k_{tr}\).

The dosing records will be associated with the first transit compartment, in this case, CMT1. This compartment does not count towards \(N_{cpt}\). The distinct absorption compartment which deposits into the central compartment also does not count towards \(N_{cpt}\). The PK records will be associated with the CENTRAL compartment, in this case, CMT2.

CMTs 3, 4, 5, and 6, represent the transit compartment chain (4 compartments total). CMT1 will feed into CMT3, which feeds to 4, which feeds to 5, which feeds to 6. CMT6 feed into CMT7, the distinct absorption compartment. Drug is absorbed into the central compartment from the distinct absorption compartment with first order rate constant \(ka\).

;4 Transit Compartments into Distinct Absorption Compartment
;First-Order Elimination (ke) from Central
;Transit Compartment 0=CMT1, CENTRAL=CMT2
;Transit Compartments 1, 2, 3, 4 = CMT3/4/5/6 
;Distinct absorption compartment = CMT7

$INPUT IGNORE=C ID TIME DV AMT EVID CMT RATE

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DOSING)
COMP=(CENTRAL)
COMP=(TRANS1)
COMP=(TRANS2)
COMP=(TRANS3)
COMP=(TRANS4)
COMP=(ABSORB)
 
$PK
NCPT=4 ;4 Transit compartments
TVMTT = THETA(1) ;Mean Transit Time
TVKA = THETA(2) ;Absorption rate constant from the distinct absorption compartment

;Typical Values
TVVC = THETA(3)
TVCL = THETA(4)

;Parameters
MTT = TVMTT*EXP(ETA(1))
KA = TVKA*EXP(ETA(2))
VC = TVVC*EXP(ETA(3))
CL = TVCL*EXP(ETA(4))

;Derived Parameters
KTR = (NCPT+1)/MTT ;Transit rate through the linked compartments
K20 = CL/VC
K72 = KA

$DES
DADT(1) = -KTR*A(1) ;Dosing compartment
DADT(3) = KTR*(A(1) - A(3)) ;Transit compartment 1
DADT(4) = KTR*(A(3) - A(4)) ;Transit compartment 2
DADT(5) = KTR*(A(4) - A(5)) ;Transit compartment 3
DADT(6) = KTR*(A(5) - A(6)) ;Transit compartment 4
DADT(7) = KTR*A(6) - K72*A(7) ;Distinct absorption compartment
DADT(2) = K72*A(7) - K20*A(2) ;Central PK compartment

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Stirling Approximation to Estimate the Number of Transit Compartments

This transit model cleverly utilizes a continuous function that closely resembles the factorial function (Stirling Approximation), shown below.

\(N! = \sqrt{2\pi}*N^{N+0.5}*e^{-N}\)

The advantage of this approach is that \(N_{cpt}\), the number of transit compartments, is no longer restricted to integer values and can now be estimated using maximum likelihood methods, eliminating the iterative stepwise testing of \(N_{cpt}\) associated with the Erlang distribution model. The analyst may also elect to add inter-individual variability on \(N_{cpt}\).

This model can also be used as a launching pad to quickly determine the ideal number of transit compartments for explicit encoding in a separate model. This control stream will work for multiple doses, with the assumption that virtually all of the drug has been absorbed into the central compartment prior to the next dose.

For linked transit compartments, \(k_{tr} = \frac{N_{cpt} + 1}{MTT}\), where \(N_{cpt}\) is the number of transit compartments, not including the initial compartment which initializes the bolus dose.

Recall the analytical solution for the instantaneous rate of compound exiting the final transit compartment for integer values of \(N_{cpt}\) is \(f(t) = Dose*\frac{k_{tr}^{N_{cpt} + 1}}{\Gamma{}(N_{cpt} + 1)}*t^{N_{cpt}}*e^{-k_{tr}t}\)

Given that \(\Gamma{}(N_{cpt} + 1) = N_{cpt}!\), we can apply the Stirling approximation to obtain an approximate solution for \(f(t)\)

\(f(t) = Dose*\frac{k_{tr}^{N_{cpt} + 1}}{\sqrt{2\pi}~*~N_{cpt}^{(N_{cpt}+0.5)}~*~e^{-N_{cpt}}}*t^{N_{cpt}}*e^{-k_{tr}t}\)

In this control stream, the final transit compartment directly deposits into the central compartment. Since the solution is analytical, the value of CMT1 in the \(\$DES\) block does not matter.

The dosing records will be associated with the DEPOT compartment, in this case, CMT1. The PK records will be associated with the CENTRAL compartment, in this case, CMT2.

;Stirling Approximation to estimate number of transit compartments
;First-Order Elimination (ke) from Central
;DEPOT=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT RATE

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Update Doses and Dose Times
IF(AMT.GT.0)DOSE = AMT
IF(AMT.GT.0)TDOSE = TIME

;Erlang Inputs
TVNCPT = THETA(1) ;Estimate number of compartments
TVMTT = THETA(2) ;Mean Transit Time
F1 = 0; We do not use the AMT in the dosing record as the input is analytically derived

;Typical Values
TVVC = THETA(3)
TVCL = THETA(4)

;Parameters

NCPT = TVNCPT*EXP(ETA(1))
MTT = TVMTT*EXP(ETA(2))
VC = TVVC*EXP(ETA(3))
CL = TVCL*EXP(ETA(4))

;Derived Parameters
KTR = (NCPT+1)/MTT ;Transit rate through the linked compartments
K20 = CL/VC

;Approximation of NCPT factorial (use if using identity scale expression for DADT(2))
;SAPPX = SQRT(2*3.1415)*NCPT**(NCPT+0.5)*EXP(-NCPT) 

;Log-transformation of NCPT factorial (use if using log transform expression for DADT(2))
LNSAPPX = LOG(2.5066) + (NCPT+0.5)*LOG(NCPT) - NCPT 


$DES
TCUR = T - TDOSE ;Time after dosing
DADT(1) = 0
;DADT(2) = DOSE*(KTR**(NCPT+1))*(TCUR**NCPT)*EXP(-KTR*TCUR)/SAPPX  - K20*A(2) ;Identity Scale (if this causes numerical difficulty use the log-transform)
DADT(2) = EXP(PLOG(DOSE) + (NCPT+1)*PLOG(KTR) + NCPT*PLOG(TCUR) - KTR*TCUR - LNSAPPX) - K20*A(2) ;Log-transform expression, use PLOG to protect against LOG(0)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Stirling Approximation (Into a Distinct Absorption Compartment)

This model can also be used as a launching pad to quickly determine the ideal number of transit compartments for explicit encoding in a separate model. This control stream will work for multiple doses, with the assumption that virtually all of the drug has been absorbed into the central compartment prior to the next dose.

;Stirling Approximation to estimate number of transit compartments (into a distinct absorption compartment)
;First-Order Elimination (ke) from Central
;DEPOT=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT RATE

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Update Doses and Dose Times
IF(AMT.GT.0)DOSE = AMT
IF(AMT.GT.0)TDOSE = TIME

;Erlang Inputs
TVNCPT = THETA(1) ;Estimate number of compartments
TVMTT = THETA(2) ;Mean Transit Time
TVKA = THETA(3) ; absorption rate constant from the DEPOT/CMT1 into CENTRAL (distinct from ktr)
F1 = 0; We do not use the AMT in the dosing record as the input into CMT1 is defined analytically

;Typical Values
TVVC = THETA(4)
TVCL = THETA(5)

;Parameters

NCPT = TVNCPT*EXP(ETA(1))
MTT = TVMTT*EXP(ETA(2))
KA = TVKA*EXP(ETA(3))
VC = TVVC*EXP(ETA(4))
CL = TVCL*EXP(ETA(5))

;Derived Parameters
KTR = (NCPT+1)/MTT ;Transit rate through the linked compartments
K12 = KA ;absorption rate constant from DEPOT to CENTRAL
K20 = CL/VC

;Approximation of NCPT factorial (use if using identity scale expression for DADT(2))
;SAPPX = SQRT(2*3.1415)*NCPT**(NCPT+0.5)*EXP(-NCPT) 

;Log-transformation of NCPT factorial (use if using log transform expression for DADT(2))
LNSAPPX = LOG(2.5066) + (NCPT+0.5)*LOG(NCPT) - NCPT 


$DES
TCUR = T - TDOSE ;Time after dosing

;DADT(1) = DOSE*(KTR**(NCPT+1))*(TCUR**NCPT)*EXP(-KTR*TCUR)/SAPPX  - K12*A(1) ;Identity Scale (if this causes numerical difficulty use the log-transform)
DADT(1) = EXP(PLOG(DOSE) + (NCPT+1)*PLOG(KTR) + NCPT*PLOG(TCUR) - KTR*TCUR - LNSAPPX) - K12*A(1) ;Log-transform expression, use PLOG to protect against LOG(0)
DADT(2) = K12*A(1) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Empirical Models (Weibull, Inverse-Gamma, Time-Varying Ka)

The mixed zero-order and first-order models above can frequently capture unusual absorption profiles and have the added benefit of intuitively associating with human physiology. However, sometimes the absorption profile is best captured empirically through more complicated but highly flexible mathematical functions, or with absorption rates that change with time.

The Weibull function is a go-to when describing in vitro dissolution profiles of compounds. It has been adapted for pharmacokinetic purposes.

Single Weibull

This control stream will also work for multiple doses, with the assumption that virtually all of the drug has been absorbed from the depot prior to the next dose.

Notice that the amount in the DEPOT A(1) is never used in \(\$DES\). The rate of transfer between compartments is only dependent on dose, time since the last dose, and the shape parameters alpha and beta. It is possible for single-dose modeling to write the \(\$DES\) block as a function of A(1). But the below control stream works for both single and multiple doses.

;Single Weibull function Absorption
;First-Order Elimination (ke) from Central
;DEPOT=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Update Doses and Dose Times
IF(AMT.GT.0)DOSE = AMT
IF(AMT.GT.0)TDOSE = TIME

F1 = 0; The AMT is not used as the input into the central compartment is analytically defined

;Typical Values
TVALPHA = THETA(1)
TVBETA = THETA(2)
TVVC = THETA(3)
TVCL = THETA(4)

;Parameters
ALPHA = TVALPHA*EXP(ETA(1))
BETA = TVBETA*EXP(ETA(2))
VC = TVVC*EXP(ETA(3))
CL = TVCL*EXP(ETA(4))

;Derived Parameters
K20 = CL/VC

$DES
TCUR = T - TDOSE

DADT(1) = 0
DADT(2) = DOSE*(BETA/ALPHA)*(TCUR/ALPHA)**(BETA - 1)*EXP(-(TCUR/ALPHA)**BETA) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Double Weibull

The same assumptions apply for this control stream as in the single Weibull control stream. IIV on WB1 should include a logit-transformation. No lag times are needed. Probably overparamaterized.

;Double Weibull function Absorption
;First-Order Elimination (ke) from Central
;DEPOT=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Update Doses and Dose Times
IF(AMT.GT.0)DOSE = AMT
IF(AMT.GT.0)TDOSE = TIME

F1 = 0;  The AMT is not used as the input into the central compartment is analytically defined

;Typical Values
TVALPHA1 = THETA(1)
TVBETA1 = THETA(2)
TVALPHA2 = THETA(3)
TVBETA2 = THETA(4)
TVWB1 = THETA(5) ;Fraction delivered via Weibull Function 1
TVVC = THETA(6)
TVCL = THETA(7)

;Parameters
ALPHA1 = TVALPHA1*EXP(ETA(1))
BETA1 = TVBETA1*EXP(ETA(2))
ALPHA2 = TVALPHA2*EXP(ETA(3))
BETA2 = TVBETA2*EXP(ETA(4))
WB1 = TVWB1 ;Fraction delivered via Weibull Function 1

VC = TVVC*EXP(ETA(5))
CL = TVCL*EXP(ETA(6))

WB2 = 1 - WB1 ;Fraction delivered via Weibull Function 2

;Derived Parameters
K20 = CL/VC

$DES
TCUR = T - TDOSE
WBRATE1 = WB1*DOSE*(BETA1/ALPHA1)*(TCUR/ALPHA1)**(BETA1 - 1)*EXP(-(TCUR/ALPHA1)**BETA1)
WBRATE2 = WB2*DOSE*(BETA2/ALPHA2)*(TCUR/ALPHA2)**(BETA2 - 1)*EXP(-(TCUR/ALPHA2)**BETA2)

DADT(1) = 0
DADT(2) =  WBRATE1 + WBRATE2 - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

The inverse gamma density function is not commonly used but is mentioned in review articles about absorption modeling in pharmacometrics.

Inverse Gamma Density Function

This control stream will also work for multiple doses, with the assumption that virtually all of the drug has been absorbed from the depot prior to the next dose.

Notice that the amount in the DEPOT A(1) is never used in \(\$DES\). The rate of transfer between compartments is only dependent on dose, time since the last dose, MAT (mean absorption time), and NVSQ (normalized variance of the gaussian density function).

;Inverse Gamma Density Absorption Function
;First-Order Elimination (ke) from Central
;DEPOT=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Update Doses and Dose Times
IF(AMT.GT.0)DOSE = AMT
IF(AMT.GT.0)TDOSE = TIME

F1 = 0;  The AMT is not used as the input into the central compartment is analytically defined

;Typical Values
TVMAT = THETA(1)
TVNVSQ = THETA(2)
TVVC = THETA(3)
TVCL = THETA(4)

;Parameters
MAT = TVMAT*EXP(ETA(1))
NVSQ = TVNVSQ*EXP(ETA(2))
VC = TVVC*EXP(ETA(3))
CL = TVCL*EXP(ETA(4))

;Derived Parameters
K20 = CL/VC

$DES
TCUR = T - TDOSE
A = SQRT(MAT/(2*3.14159*NVSQ*(TCUR**3)) ;to simplify inv gamma expression
B = ((TCUR - MAT)**2)/(2*NVSQ*MAT*TCUR) ;to simplify inv gamma expression
RATEINVGAM = DOSE*A*EXP(-B)

DADT(1) = 0
DADT(2) =  RATEINVGAM - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

The control streams below reflect models where the absorption rate begins at a baseline value and then increases or decreases to another value as a function of time. This could reflect some induction or suppression of absorption over time. This type of time-dependent change in a parameter value is more frequently encountered on clearance parameters.

Ka Sigmoidal, Changes with Time

In this control stream, the ka does not revert to the original value upon each dose. The Hill factor, GAMMA, can be fixed to 1 in the \(\$THETA\) block to remove sigmoidal character from the time course of ka (resulting in a hyperbolic ka profile).

KABASE is the baseline absorption rate at dosing, and KAIND is the added induced absorption rate which appears after dosing has occured. The asymptotic absorption rate reached is the sum of KABASE and KAIND. If KAIND < 0, then the absorption rate is decreasing with time, and if KAIND > 0, the absorption rate is increasing with time.

\(k_{a}(t) = k_{a,base} + k_{a,induced}*\frac{t^{\gamma{}}}{T^{\gamma{}}_{50} + t^\gamma{}}\)

;Sigmoidal Absorption Rate KA changes with time
;First-Order Elimination (ke) from Central
;DEPOT=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK

;Typical Values
TVKABASE = THETA(1) ;Initial KA value
TVKAIND = THETA(2) ;KA induced over time
TVKAT50 = THETA(3) ;Time at which the KA is halfway between the intial and final value
TVGAMMA = THETA(4) ;Hill factor to influence sigmoidicity
TVVC = THETA(5)
TVCL = THETA(6)

;Parameters
KABASE = TVKABASE*EXP(ETA(1))
KAIND = TVKAIND + ETA(2) ;Not necessarily positive though be cautious of negative absorptions
KAT50 = TVKAT50
GAMMA = TVGAMMA
VC = TVVC*EXP(ETA(3))
CL = TVCL*EXP(ETA(4))

;Derived Parameters
K20 = CL/VC
KASS = KABASE + KAIND ;Steady-state asymptotic KA

$DES
KA = KABASE + KAIND*(T**GAMMA)/(KAT50**GAMMA + T**GAMMA)

DADT(1) = -KA*A(1)
DADT(2) =  KA*A(1) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Ka Sigmoidal, Decreases or Increases with Time by Subject

In this control stream, the ka does not revert to the original value upon each dose. The Hill factor, GAMMA, can be fixed to 1 in the \(\$THETA\) block to remove sigmoidal character from the time course of ka (resulting in a hyperbolic ka profile).

By exponentiating the sigmoidal expression, the CL at time 0 can be set to the baseline CL, and interindividual variability can be normally distributed on Emax such that negative Emax values correspond to inhibition of ka over time and positive Emax values correspond to augmentation of ka over time in the same patient population. In this case, an individuals \(e^{E_{max}}\) reflects the proportional change of ka from the baseline over time (i.e., if \(e^{E_{max,i}} = 2\), then subject i’s ka will asymptotically reach a value double their baseline). See Bajaj et al. (2017) CPT: Pharmacometrics and Systems Pharmacology.

\(k_{a}(t) = k_{a,base}*exp(\frac{E_{max}~*~{t^{\gamma{}}}}{T^{\gamma{}}_{50}~+~t^\gamma{}})\)

;Sigmoidal Absorption Rate, Exponentiatied  
;First-Order Elimination (ke) from Central
;DEPOT=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK

;Typical Values
TVKABASE = THETA(1) ;Initial KA value
TVEMAX = THETA(2) ;EMAX/IMAX for sigmoidal absorption time profile
TVKAT50 = THETA(3) ;Time at which the KA is halfway between the intial and final value
TVGAMMA = THETA(4) ;Hill factor to influence sigmoidicity
TVVC = THETA(5)
TVCL = THETA(6)

;Parameters
KABASE = TVKABASE*EXP(ETA(1))
EMAX = TVEMAX + ETA(2) ;Normally distributed IIV
KAT50 = TVKAT50
GAMMA = TVGAMMA
VC = TVVC*EXP(ETA(3))
CL = TVCL*EXP(ETA(4))

;Derived Parameters
K20 = CL/VC

$DES
KA = KABASE*EXP((EMAX*(T**GAMMA))/(KAT50**GAMMA + T**GAMMA))

DADT(1) = -KA*A(1)
DADT(2) =  KA*A(1) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Ka Changes with Time Exponentially

In this control stream, the ka does not revert to the original value upon each dose. KABASE is the baseline absorption rate at dosing, and KAIND is the added induced absorption rate after dosing has occured. The asymptotic absorption rate reached is the sum of KABASE and KAIND. If KAIND < 0, then the absorption rate is decreasing with time, and if KAIND > 0, the absorption rate is increasing with time.

\(k_{a}(t) = k_{a,base} + k_{a,induced}*(1 - e^{-IND~*~t})\)

;Absorption Rate Changes with Time (Exponentially)
;First-Order Elimination (ke) from Central
;DEPOT=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK

;Typical Values
TVKABASE = THETA(1) ;Initial KA value
TVKAIND = THETA(2) ;Induced KA value
TVIND = THETA(3) ;Induction rate constant
TVVC = THETA(4)
TVCL = THETA(5)

;Parameters
KABASE = TVKABASE*EXP(ETA(1))
KAIND = TVKAIND + ETA(2) ;Not necessarily positive though be cautious of negative absorptions
IND = TVIND
VC = TVVC*EXP(ETA(3))
CL = TVCL*EXP(ETA(4))

;Derived Parameters
K20 = CL/VC

$DES
KA = KABASE + KAIND*(1 - EXP(-IND*T))

DADT(1) = -KA*A(1)
DADT(2) =  KA*A(1) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Ka Changes with Time (Weibull-Like)

In this control stream, the ka does not revert to the original value upon each dose. KABASE is the baseline absorption rate at dosing, and KAIND is the added induced absorption rate after dosing has occured. The asymptotic absorption rate reached is the sum of KABASE and KAIND. If KAIND < 0, then the absorption rate is decreasing with time, and if KAIND > 0, the absorption rate is increasing with time.

\(k_{a}(t) = k_{a,base} + k_{a,induced}*(1 - e^{-(k~*~t)^{\gamma{}}})\)

;Absorption Rate Changes with Time (Exponentially)
;First-Order Elimination (ke) from Central
;DEPOT=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK

;Typical Values
TVKABASE = THETA(1) ;Initial KA value
TVKAIND = THETA(2) ;Induced KA value
TVWBK = THETA(3) ;Weibull parameter K
TVWBGAMMA = THETA(4) ;Weibull parameter Gamma
TVVC = THETA(5)
TVCL = THETA(6)

;Parameters
KABASE = TVKABASE*EXP(ETA(1))
KAIND = TVKAIND + ETA(2) ;Not necessarily positive though be cautious of negative absorptions
WBK = TVWBK
WBGAMMA = TVWBGAMMA 
VC = TVVC*EXP(ETA(3))
CL = TVCL*EXP(ETA(4))

;Derived Parameters
K20 = CL/VC

$DES
KA = KABASE + KAIND*(1 - EXP(-(WBK*T)**WBGAMMA))

DADT(1) = -KA*A(1)
DADT(2) =  KA*A(1) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Reverting Ka to Baseline At Each Dose

In the above time-varying ka control streams, the ka value does not return to its baseline value after each dose. Once changed, it remains changed. To have the ka value reset to its initial value at each dosing occasion, the following changes can be made to the control stream. Ideally, the DEPOT compartment will be empty or close to empty between dosing occasions.

First, using IF statements, save the time of each dose within \(\$PK\). Then, within \(\$DES\), define the current time since the most recent dose TCUR.

Finally, use TCUR for T when T appears in the time-dependent expressions for KA. Do not a pre-made column TAD in the dataset. TCUR should be calculated in \(\$DES\) to ensure fine-grained sampling of times for continuous calculation of \(ka(t)\).

$PK
;Update Doses and Dose Times
IF(AMT.GT.0)TDOSE = TIME

$DES
TCUR = T - TDOSE
;Then replace T with TCUR in ka = f(T)

Distribution Models

For pharmacokinetic modeling, selecting between options for modeling the distribution of drug in the body generally refers to the number of theoretical compartments that the drug can distribute between. Conceptually, the drug in plasma can equilibriate with other bodily tissues like fat. For each additional theoretical compartment that drug can transfer into from plasma, an intercompartmental clearance \(Q\) and peripheral volume \(V_p\) is declared.

Let the central plasma PK compartment be CMT2 with volume \(V_c\) and the peripheral compartment be CMT3 with volume \(V_p\).

The mass transfer from CMT2 into CMT3 can be defined as \(k_{23}*A(2)\) where \(k_{23}\) is the transfer rate constant and \(A(2)\) is the amount of drug in CMT2. \(k_{23}\) can be parameterized as \(\frac{Q}{V_c}\).

The mass transfer in the reverse direction, from CMT3 into CMT2, can be defined as \(k_{32}*A(3)\). \(k_{32}\) can be parameterized as \(\frac{Q}{V_p}\).

For simplicity, the one, two, and three compartment models shown below are first-order absorption models which have first-order elimination from the central compartment.

One-Compartment Model

The dosing records will be associated with the DEPOT compartment, in this case, CMT1. The PK records will be associated with the CENTRAL compartment, in this case, CMT2.

$PROBLEM One-Compartment Model
;First-Order Absorption (ka) from Depot into Central
;First-Order Elimination (ke) from Central
;DEPOT=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT 

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Typical Values
TVKA = THETA(1)
TVVC = THETA(2)
TVCL = THETA(3)

;Parameters
KA = TVKA*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
CL = TVCL*EXP(ETA(3))

;Derived Parameters
K12 = KA
K20 = CL/VC

$DES
DADT(1) = -K12*A(1)
DADT(2) =  K12*A(1) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Two-Compartment Model

The dosing records will be associated with the DEPOT compartment, in this case, CMT1. The PK records will be associated with the CENTRAL compartment, in this case, CMT2. CMT3 is the theoretical peripheral compartment.

$PROBLEM Two-Compartment Model
;First-Order Absorption (ka) from Depot into Central
;Distribution between Central and Peripheral Volume
;First-Order Elimination (ke) from Central
;DEPOT=CMT1, CENTRAL=CMT2, PERIPH=CMT3

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
COMP=(PERIPH)
 
$PK
;Typical Values
TVKA = THETA(1)
TVVC = THETA(2)
TVCL = THETA(3)
TVVP = THETA(4)
TVQ = THETA(5)

;Parameters
KA = TVKA*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
CL = TVCL*EXP(ETA(3))
VP = TVVP
Q = TVQ

;Derived Parameters
K12 = KA
K23 = Q/VC
K32 = Q/VP
K20 = CL/VC

$DES
DADT(1) = -K12*A(1)
DADT(2) =  K12*A(1) + K32*A(3) - K23*A(2) - K20*A(2)
DADT(3) =  K23*A(2) - K32*A(3)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Three-Compartment Model

The dosing records will be associated with the DEPOT compartment, in this case, CMT1. The PK records will be associated with the CENTRAL compartment, in this case, CMT2. CMT3 is the theoretical peripheral compartment 1. CMT3 is the theoretical peripheral compartment 2.

$PROBLEM Three-Compartment Model
;First-Order Absorption (ka) from Depot into Central
;Distribution between Central and Two Peripheral Volumes (VP1 and VP2)
;First-Order Elimination (ke) from Central
;DEPOT=CMT1, CENTRAL=CMT2, PERIPH=CMT3

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
COMP=(PERIPH)
 
$PK
;Typical Values
TVKA = THETA(1)
TVVC = THETA(2)
TVCL = THETA(3)
TVVP1 = THETA(4)
TVQ1 = THETA(5)
TVVP2 = THETA(6)
TVQ2 = THETA(7)

;Parameters
KA = TVKA*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
CL = TVCL*EXP(ETA(3))
VP1 = TVVP1
Q1 = TVQ1
VP2 = TVVP2
Q2 = TVQ2

;Derived Parameters
K12 = KA
K23 = Q1/VC
K32 = Q1/VP1
K24 = Q2/VC
K42 = Q2/VP2
K20 = CL/VC

$DES
DADT(1) = -K12*A(1)
DADT(2) =  K12*A(1) + K32*A(3) + K42*A(4) - K23*A(2) - K24*A(2) - K20*A(2)
DADT(3) =  K23*A(2) - K32*A(3)
DADT(4) =  K24*A(2) - K42*A(4)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Elimination Models

For simpler pharmacokinetic modeling (without pharmacodynamic considerations), the options for describing the elimination of drug from the plasma compartment are generally first-order or saturable (Michaelis-Menten) clearance mechanisms.

The kidneys and liver can filter out drug from a certain volume of plasma per unit time. This lends itself nicely to modeling the rate of drug removal from plasma as proportional to the concentration of drug in the plasma (i.e. first-order elimination). In cases where the drug must be removed from the plasma through binding of a receptor (perhaps on the surface of hepatocytes), or if the drug is metabolized by enzymes that can be saturated at doses used in the study, Michaelis-Menten elimination can be applied.

Linear and Saturable Clearance

First-Order Elimination

The dosing records will be associated with the DEPOT compartment, in this case, CMT1. The PK records will be associated with the CENTRAL compartment, in this case, CMT2.

At the TIME of dosing records, an instantaneous amount AMT will be added to the DEPOT compartment.

;First-Order Absorption (ka) from Depot into Central
;First-Order Elimination (k20) from Central
;DEPOT=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Typical Values
TVKA = THETA(1)
TVVC = THETA(2)
TVCL = THETA(3)

;Parameters
KA = TVKA*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
CL = TVCL*EXP(ETA(3))

;Derived Parameters
K12 = KA
K20 = CL/VC

$DES
DADT(1) = -K12*A(1)
DADT(2) =  K12*A(1) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Michaelis-Menten Elimination

The dosing records will be associated with the DEPOT compartment, in this case, CMT1. The PK records will be associated with the CENTRAL compartment, in this case, CMT2.

At the TIME of dosing records, an instantaneous amount AMT will be added to the DEPOT compartment.

;First-Order Absorption (ka) from Depot into Central
;Saturable (Michaelis-Menten Elimination) (Vmax/KM) from Central
;DEPOT=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Typical Values
TVKA = THETA(1)
TVVC = THETA(2)
TVVMAX = THETA(3)
TVKM = THETA(4)

;Parameters
KA = TVKA*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
VMAX = TVVMAX*EXP(ETA(3))
KM = TVKM*EXP(ETA(4))

;Derived Parameters
K12 = KA

$DES
C2 = A(2)/VC ;Be cautious with units to have them match the desired units of KM

DADT(1) = -K12*A(1)
DADT(2) =  K12*A(1) - (VMAX*C2)/(KM + C2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Additionally, first-order and Michaelis-Menten elimination can simultaneously be present if drug is removed via two these two pathways in parallel. However, the data collection during the elimination phase will likely need to be intensively sampled to accurately estimate parameters for dual clearance.

Dual First-Order and Michaelis-Menten Elimination

The dosing records will be associated with the DEPOT compartment, in this case, CMT1. The PK records will be associated with the CENTRAL compartment, in this case, CMT2.

At the TIME of dosing records, an instantaneous amount AMT will be added to the DEPOT compartment.

;First-Order Absorption (ka) from Depot into Central
;Dual First-Order (k20) and Saturable (Michaelis-Menten Elimination) (Vmax/KM) from Central
;DEPOT=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT RATE

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Typical Values
TVKA = THETA(1)
TVVC = THETA(2)
TVCL = THETA(3)
TVVMAX = THETA(4)
TVKM = THETA(5)

;Parameters
KA = TVKA*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
CL = TVCL*EXP(ETA(3))
VMAX = TVVMAX*EXP(ETA(4))
KM = TVKM*EXP(ETA(5))

;Derived Parameters
K12 = KA
K20 = CL/VC

$DES
C2 = A(2)/VC ;Be cautious with units to have them match the desired units of KM

DADT(1) = -K12*A(1)
DADT(2) =  K12*A(1) - K20*A(2) - (VMAX*C2)/(KM + C2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Time-Varying Clearance

In some cases, especially with drugs that may induce or repress the expression of enzymes involved in drug metabolism (CYP family), the clearance of drug can vary with time. Below are a variety of control streams where the clearance is encoded as a function of time. Their structure is very similar to the time-varying absorption rates described in an earlier section.

CL Sigmoidal, Changes with Time

In this control stream, the CL does not revert to the original value upon each dose. The Hill factor, GAMMA, can be fixed to 1 in the \(\$THETA\) block to remove sigmoidal character from the time course of CL (resulting in a hyperbolic CL profile).

CLBASE is the baseline clearance at dosing, and CLIND is the added clearance after dosing has occurred. The asymptotic clearance reached is the sum of CLBASE and CLIND. If CLIND < 0, then the clearance decreases with time, and if CLIND > 0, the clearance is increasing with time (more common).

\(CL(t) = CL_{base} + CL_{induced}*\frac{t^{\gamma{}}}{T^{\gamma{}}_{50} + t^\gamma{}}\)

;Sigmoidal clearance CL changes with time
;First-Order Elimination (ke) from Central
;DEPOT=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT 

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK

;Typical Values
TVKA = THETA(1)
TVVC = THETA(2)
TVCLBASE = THETA(3) ;Initial CL value
TVCLIND = THETA(4) ;CL induced over time
TVCLT50 = THETA(5) ;Time at which the CL is halfway between the intial and final value
TVGAMMA = THETA(6) ;Hill factor to influence sigmoidicity

;Parameters
KA = TVKA*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
CLBASE = TVCLBASE*EXP(ETA(3))
CLIND = TVCLIND + ETA(4) ;Not necessarily positive though be cautious of negative clearances
CLT50 = TVCLT50
GAMMA = TVGAMMA

;Derived Parameters
CLSS = CLBASE + CLIND ;Steady-state asymptotic CL

$DES
CL = CLBASE + CLIND*(T**GAMMA)/(CLT50**GAMMA + T**GAMMA)
K20 = CL/VC

DADT(1) = -KA*A(1)
DADT(2) =  KA*A(1) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

CL Sigmoidal, Decreases or Increases with Time by Subject

In this control stream, the CL does not revert to the original value upon each dose. The Hill factor, GAMMA, can be fixed to 1 in the \(\$THETA\) block to remove sigmoidal character from the time course of CL (resulting in a hyperbolic CL profile).

By exponentiating the sigmoidal expression, the CL at time 0 can be set to the baseline CL, and interindividual variability can be normally distributed on Emax such that negative Emax values correspond to inhibition of CL over time and positive Emax values correspond to augmentation of CL over time in the same patient population. In this case, an individuals \(e^{E_{max}}\) reflects the proportional change of CL from the baseline over time (i.e., if \(e^{E_{max,i}} = 2\), then subject i’s CL will asymptotically reach a value double their baseline). See Bajaj et al. (2017) CPT: Pharmacometrics and Systems Pharmacology.

\(CL(t) = CL_{base}*exp(\frac{E_{max}~*~{t^{\gamma{}}}}{T^{\gamma{}}_{50}~+~t^\gamma{}})\)

;Sigmoidal CL, Exponentiatied  
;First-Order Elimination (ke) from Central
;DEPOT=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK

;Typical Values
TVKA = THETA(1)
TVVC = THETA(2)
TVCLBASE = THETA(3) ;Initial CL value
TVEMAX = THETA(4) ;EMAX/IMAX for sigmoidal CL time profile
TVCLT50 = THETA(5) ;Time at which the CL is halfway between the intial and final value
TVGAMMA = THETA(6) ;Hill factor to influence sigmoidicity

;Parameters
KA = TVKA*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
CLBASE = TVCLBASE*EXP(ETA(3))
EMAX = TVEMAX + ETA(4) ;Normally distributed IIV
CLT50 = TVCLT50
GAMMA = TVGAMMA

$DES
CL = CLBASE*EXP((EMAX*(T**GAMMA))/(CLT50**GAMMA + T**GAMMA))
K20 = CL/VC

DADT(1) = -KA*A(1)
DADT(2) =  KA*A(1) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

CL Changes with Time Exponentially

In this control stream, the CL does not revert to the original value upon each dose. CLBASE is the baseline clearance at dosing, and CLIND is the added induced clearance after dosing has occured. The asymptotic clearance reached is the sum of CLBASE and CLIND. If CLIND < 0, then the clearance is decreasing with time, and if CLIND > 0, the clearance is increasing with time.

\(CL(t) = CL_{base} + CL_{induced}*(1 - e^{-IND~*~t})\)

;Clearance Changes with Time (Exponentially)
;First-Order Elimination (ke) from Central
;DEPOT=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT RATE

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK

;Typical Values
TVKA = THETA(1)
TVVC = THETA(2) 
TVCLBASE = THETA(3) ;Baseline CL
TVCLIND = THETA(4) ;Induced CL
TVIND = THETA(5) ;CL Induction Rate Constant

;Parameters
KA = TVKABASE*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
CLBASE = TVCLBASE*EXP(ETA(3))
CLIND = TVCLIND + ETA(4) ;Not necessarily positive but be cautious with negative clearances
IND = TVIND

;Derived Parameters
CLSS = CLBASE + CLIND ;Asymptotic CL 

$DES
CL = CLBASE + CLIND*(1 - EXP(-IND*T))
K20 = CL/VC

DADT(1) = -KA*A(1)
DADT(2) =  KA*A(1) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

CL Changes with Time (Weibull-Like)

In this control stream, the CL does not revert to the original value upon each dose. CLBASE is the baseline clearance at dosing, and CLIND is the added induced clearance after dosing has occured. The asymptotic clearance reached is the sum of CLBASE and CLIND. If CLIND < 0, then the clearance is decreasing with time, and if CLIND > 0, the clearance is increasing with time.

\(CL(t) = CL_{base} + CL_{induced}*(1 - e^{{-(k~*~t)}^{\gamma{}}})\)

;Clearance Changes with Time (Weibull-Like)
;First-Order Elimination (ke) from Central
;DEPOT=CMT1, CENTRAL=CMT2

$INPUT IGNORE=C ID TIME DV AMT EVID CMT RATE

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK

;Typical Values
TVKA = THETA(1)
TVVC = THETA(2) 
TVCLBASE = THETA(3) ;Baseline CL
TVCLIND = THETA(4) ;Induced CL
TVWBK = THETA(5) ;Weibull K Parameter
TVWBGAMMA = THETA(6) ;Weibull Gamma Parameter

;Parameters
KA = TVKABASE*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
CLBASE = TVCLBASE*EXP(ETA(3))
CLIND = TVCLIND + ETA(4) ;Not necessarily positive but be cautious with negative clearances
WBK = TVWBK
WBGAMMA = TVWBGAMMA

;Derived Parameters
CLSS = CLBASE + CLIND ;Asymptotic CL 

$DES
CL = CLBASE + CLIND*(1 - EXP((-WBK*T)**WBGAMMA))
K20 = CL/VC

DADT(1) = -KA*A(1)
DADT(2) =  KA*A(1) - K20*A(2)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Reverting CL to Baseline At Each Dose

In the above time-varying CL control streams, the CL value does not return to its baseline value after each dose. Once changed, it remains changed. To have the CL value reset to its initial value at each dosing occasion, the following changes can be made to the control stream. Ideally, most of the drug will have been cleared between doses.

First, using IF statements, save the time of each dose within \(\$PK\). Then, within \(\$DES\), define current time since the most recent dose TCUR.

Finally, use TCUR for \(T\) when \(T\) appears in the time-dependent expressions for CL. Do not a pre-made column TAD in the dataset. TCUR should be calculated in \(\$DES\) to ensure fine-grained sampling of times for continuous calculation of \(CL(t)\).

$PK
;Update Doses and Dose Times
IF(AMT.GT.0)TDOSE = TIME

$DES
TCUR = T - TDOSE
;Replace T with TCUR for all appearances in CL = f(T)

Enzyme Induction Model

Perhaps the drug induces expression of enzymes that facilitate clearance with a significant delay between the plasma concentration of drug and the enzyme increase in concentration. In this case, we may elect to use an virtual enzyme turnover model. There are many more parameters to estimate in this model than when the clearance is expressed as a function of time. In the case of induction, the plasma concentration of drug will augment the zero-order synthesis rate of “enzyme” in a virtual effect compartment with the steady-state initial “enzyme concentration” of 1. After some time, the effect will increase based on the presence of drug in the plasma.

The clearance of drug from the central compartment will then be a function (proportional is simplest) of the instantaneous value in the effect compartment. Conceptually, this is a turnover model where the virtual enzyme compartment is the PD response (and the PD response is proportional to the clearance).

Induction

When using a linear clearance,

\(CL(t) = CL_{base}*A_{enz}\)

The amount in the enzyme effect compartment is: \(\frac{dA_{enz}}{dt} = k_{enz}*(1 + \frac{E_{max}~*~C_{p}}{EC_{50} + C_{p}}) - k_{enz}*A_{enz}\)

where \(C_{p}\) is the current plasma concentration \(\frac{A(2)}{V_c}\) and the initial condition of \(A_{enz}(t=0) = 1\).

A saturable clearance can also be used, where

\(CL(t) = \frac{V_{max}*C_{p}}{K_{m} + C_{p}}*A_{enz}\). The differential equation describing \(A_{enz}\) is the same as above.

Enzyme Induction of Clearance, Linear Clearance

;First-Order Absorption (ka) into central
;Linear Clearance from central compartment with enzyme auto-induction
;DEPOT=CMT1, CENTRAL=CMT2, ENZYME EFFECT=CMT3

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
COMP=(ENZ)
 
$PK

;Typical Values
TVKA = THETA(1)
TVVC = THETA(2)
TVCL = THETA(3) ;Baseline Clearance
TVKENZ = THETA(4) ;Zero-order synthesis rate of enzyme
TVEMAX = THETA(5) ;Maximum effect of drug on enzyme synthesis rate
TVEC50 = THETA(6) ;Half-maximal concentration of drug on enzyme synthesis rate

;Parameters
KA = TVKA*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
CL = TVCL*EXP(ETA(3))
KENZ = TVKENZ
EMAX = TVEMAX
EC50 = TVEC50

;Derived Parameters
K20 = CL/VC ;Baseline elimination rate constant from central compartment

;Initialize enzyme compartment to 1
A_0(3) = 1

$DES
C2 = A(2)/VC ;Central compartment concentration of drug

EFF = (EMAX*C2)/(EC50 + C2)

DADT(1) = -KA*A(1)
DADT(2) =  KA*A(1) - K20*A(2)*A(3)
DADT(3) =  KENZ*(1 + EFF) - KENZ*A(3)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Enzyme Induction of Clearance, Saturable Clearance

;First-Order Absorption (ka) into central
;Saturable Clearance from central compartment with enzyme auto-induction
;DEPOT=CMT1, CENTRAL=CMT2, ENZYME EFFECT=CMT3

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
COMP=(ENZ)
 
$PK

;Typical Values
TVKA = THETA(1)
TVVC = THETA(2)
TVVMAX = THETA(3) ;Maximum velocity of drug elimination
TVKM = THETA(4) ;KM of drug elimination
TVKENZ = THETA(5) ;Zero-order synthesis rate of enzyme
TVEMAX = THETA(6) ;Maximum effect of drug on enzyme synthesis rate
TVEC50 = THETA(7) ;Half-maximal concentration of drug on enzyme synthesis rate

;Parameters
KA = TVKA*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
VMAX = TVVMAX*EXP(ETA(3))
KM = TVKM*EXP(ETA(4))
KENZ = TVKENZ
EMAX = TVEMAX
EC50 = TVEC50

;Initialize enzyme compartment to 1
A_0(3) = 1

$DES
C2 = A(2)/VC ;Central compartment concentration of drug

EFF = (EMAX*C2)/(EC50 + C2)

DADT(1) = -KA*A(1)
DADT(2) =  KA*A(1) - ((VMAX*C2)/(KM + C2))*A(3)
DADT(3) =  KENZ*(1 + EFF) - KENZ*A(3)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Enzyme Inhibition Model

Perhaps the drug inhibits expression of enzymes that facilitate clearance with a significant delay between the plasma concentration of drug and the decrease in enzyme concentration. In this case, we may elect to use an virtual enzyme turnover model. There are many more parameters to estimate in this model than when the clearance is expressed as a function of time. In the case of inhibition, the plasma concentration of drug will reduce the zero-order synthesis rate of “enzyme” in a virtual effect compartment with the steady-state initial “enzyme concentration” of 1. After some time, the value in the effect compartment will be reduced based on the presence of drug in the plasma.

The clearance of drug from the central compartment will then be a function (proportional is simplest) of the instantaneous value in the effect compartment.

Inhibition

When using a linear clearance,

\(CL(t) = CL_{base}*A_{enz}\)

The amount in the enzyme effect compartment is: \(\frac{dA_{enz}}{dt} = k_{enz}*(1 - \frac{I_{max}~*~C_{p}}{IC_{50} + C_{p}}) - k_{enz}*A_{enz}\)

where \(C_{p}\) is the current plasma concentration \(\frac{A(2)}{V_c}\) and the initial condition of \(A_{enz}(t=0) = 1\).

A saturable clearance can also be used, where

\(CL(t) = \frac{V_{max}*C_{p}}{K_{m} + C_{p}}*A_{enz}\). The differential equation describing \(A_{enz}\) is the same as above.

Enzyme Inhibition of Clearance, Linear Clearance

;First-Order Absorption (ka) into central
;Linear Clearance from central compartment with enzyme auto-inhibition
;DEPOT=CMT1, CENTRAL=CMT2, ENZYME EFFECT=CMT3

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
COMP=(ENZ)
 
$PK

;Typical Values
TVKA = THETA(1)
TVVC = THETA(2)
TVCL = THETA(3) ;Baseline Clearance
TVKENZ = THETA(4) ;Zero-order synthesis rate of enzyme
TVIMAX = THETA(5) ;Maximum inhibition of enzyme synthesis rate by drug
TVIC50 = THETA(6) ;Concentration of drug that achieves half-maximal inhibition of enzyme synthesis rate

;Parameters
KA = TVKA*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
CL = TVCL*EXP(ETA(3))
KENZ = TVKENZ
IMAX = TVIMAX
IC50 = TVIC50

;Derived Parameters
K20 = CL/VC ;Baseline elimination rate constant from central compartment

;Initialize enzyme compartment to 1
A_0(3) = 1

$DES
C2 = A(2)/VC ;Central compartment concentration of drug

INH = (IMAX*C2)/(IC50 + C2)

DADT(1) = -KA*A(1)
DADT(2) =  KA*A(1) - K20*A(2)*A(3)
DADT(3) =  KENZ*(1 - EFF) - KENZ*A(3)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Enzyme Inhibition of Clearance, Saturable Clearance

;First-Order Absorption (ka) into central
;Saturable Clearance from central compartment with enzyme auto-inhibition
;DEPOT=CMT1, CENTRAL=CMT2, ENZYME EFFECT=CMT3

$INPUT IGNORE=C ID TIME DV AMT EVID CMT

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
COMP=(ENZ)
 
$PK

;Typical Values
TVKA = THETA(1)
TVVC = THETA(2)
TVVMAX = THETA(3) ;Maximum velocity of drug elimination
TVKM = THETA(4) ;KM of drug elimination
TVKENZ = THETA(5) ;Zero-order synthesis rate of enzyme
TVIMAX = THETA(6) ;Maximum inhibition of enzyme synthesis rate by drug
TVIC50 = THETA(7) ;Concentration of drug that achieves half-maximal inhibition of enzyme synthesis rate

;Parameters
KA = TVKA*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
VMAX = TVVMAX*EXP(ETA(3))
KM = TVKM*EXP(ETA(4))
KENZ = TVKENZ
IMAX = TVIMAX
IC50 = TVIC50

;Initialize enzyme compartment to 1
A_0(3) = 1

$DES
C2 = A(2)/VC ;Central compartment concentration of drug

INH = (IMAX*C2)/(IC50 + C2)

DADT(1) = -KA*A(1)
DADT(2) =  KA*A(1) - ((VMAX*C2)/(KM + C2))*A(3)
DADT(3) =  KENZ*(1 - INH) - KENZ*A(3)

;Adjust below as necessary
$ERROR
$THETA
$OMEGA
$SIGMA
$ESTIMATION
$COV 
$TABLE

Residual Error Models

Residual unexplained variability (RUV) is estimated in the \(\$\)SIGMA block. Just like \(\$\)OMEGA estimates, estimated values in \(\$\)SIGMA are considered variances and individual samples, or epsilons (sampled via EPS(N)), are taken from a normal distribution centered around 0.

However, it is common practice to FIX the \(\$\)SIGMA block to a single value of 1, and estimate RUV parameters in the \(\$\)THETA block. RUV THETA parameters are then multiplied by EPS(1) within the \(\$\)ERROR block to scale the RUV appropriately.

A useful tip to quickly compare additive, proportional, and combined additive and proportional residual error models is to encode a single combined error model and selectively FIX the additive or proportional THETAs to 0.

Additive Error Model

\(Y = F + EPS(1)\)

Proportional RUV is fixed to 0. RUV parameters ADD and PROP estimates are on the SD scale.

$PROBLEM Additive Error, First Order Absorption/Elimination

$INPUT IGNORE=C ID TIME AMT DV CMT EVID MDV

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Typical Values
TVKA = THETA(1)
TVVC = THETA(2)
TVCL = THETA(3)

;Parameters
KA = TVKA*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
CL = TVCL*EXP(ETA(3))

;Derived Parameters
K12 = KA
K20 = CL/VC

$DES
DADT(1) = -K12*A(1)
DADT(2) =  K12*A(1) - K20*A(2)

$ERROR
ADD = THETA(4)
PROP = THETA(5)

IPRED = F
W = PSQRT(ADD*ADD + PROP*PROP*IPRED*IPRED)
IF(W.LE.1E-9) W=1E-9
Y = IPRED + W*EPS(1)
IRES = DV - IPRED
IWRES = IRES/W

$THETA
(0, 1) ;KA
(0, 1) ;VC
(0, 1) ;CL
(0, 10) ;ADDERR
0 FIX ;PROPERR

$OMEGA
(0.1) ;IIV_KA
(0.1) ;IIV_VC
(0.1) ;IIV_CL

$SIGMA
1 FIX

Proportional Error Model

\(Y = F * (1 + EPS(1))\)

Additive RUV is fixed to 0. RUV parameters ADD and PROP estimates are on the SD scale.

$PROBLEM Proportional Error, First Order Absorption/Elimination

$INPUT IGNORE=C ID TIME AMT DV CMT EVID MDV

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Typical Values
TVKA = THETA(1)
TVVC = THETA(2)
TVCL = THETA(3)

;Parameters
KA = TVKA*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
CL = TVCL*EXP(ETA(3))

;Derived Parameters
K12 = KA
K20 = CL/VC

$DES
DADT(1) = -K12*A(1)
DADT(2) =  K12*A(1) - K20*A(2)

$ERROR
ADD = THETA(4)
PROP = THETA(5)

IPRED = F
W = PSQRT(ADD*ADD + PROP*PROP*IPRED*IPRED)
IF(W.LE.1E-9) W=1E-9
Y = IPRED + W*EPS(1)
IRES = DV - IPRED
IWRES = IRES/W

$THETA
(0, 1) ;KA
(0, 1) ;VC
(0, 1) ;CL
(0) FIX ;ADDERR
(0, 0.1) ;PROPERR

$OMEGA
(0.1) ;IIV_KA
(0.1) ;IIV_VC
(0.1) ;IIV_CL

$SIGMA
1 FIX

Combined Additive and Proportional Error Model

\(Y = EPS(1) + F * (1 + EPS(2))\)

Neither RUV parameter is fixed to 0. RUV parameters ADD and PROP estimates are on the SD scale.

$PROBLEM Combined Additive and Proportional Error, First Order Absorption/Elimination

$INPUT IGNORE=C ID TIME AMT DV CMT EVID MDV

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Typical Values
TVKA = THETA(1)
TVVC = THETA(2)
TVCL = THETA(3)

;Parameters
KA = TVKA*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
CL = TVCL*EXP(ETA(3))

;Derived Parameters
K12 = KA
K20 = CL/VC

$DES
DADT(1) = -K12*A(1)
DADT(2) =  K12*A(1) - K20*A(2)

$ERROR
ADD = THETA(4)
PROP = THETA(5)

IPRED = F
W = PSQRT(ADD*ADD + PROP*PROP*IPRED*IPRED)
IF(W.LE.1E-9) W=1E-9
Y = IPRED + W*EPS(1)
IRES = DV - IPRED
IWRES = IRES/W

$THETA
(0, 1) ;KA
(0, 1) ;VC
(0, 1) ;CL
(0, 0.1) ;ADDERR
(0, 0.1) ;PROPERR

$OMEGA
(0.1) ;IIV_KA
(0.1) ;IIV_VC
(0.1) ;IIV_CL

$SIGMA
1 FIX

The log-error model (often referred to as ‘log-additive’) provides convenient pharmacokinetic properties in that simulated observations cannot be negative. Like the proportional and combined additive and proportional error models, the variance of the residual error will increase with the magnitude of the prediction.

Natural-log transformed data (within the NONMEM input dataset) is used to implement this error model.

Log-Additive Error Model

\(Y = LOG(F) + EPS(1)\)

Note LNDV, natural log-transformed DV, is used as the DV column. RUV parameter LOGADD is on the SD scale.

$PROBLEM Log-Additive Error, First Order Absorption/Elimination

$INPUT IGNORE=C ID TIME AMT DV=DROP LNDV=DV CMT EVID MDV

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;Typical Values
TVKA = THETA(1)
TVVC = THETA(2)
TVCL = THETA(3)

;Parameters
KA = TVKA*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
CL = TVCL*EXP(ETA(3))

;Derived Parameters
K12 = KA
K20 = CL/VC

$DES
DADT(1) = -K12*A(1)
DADT(2) =  K12*A(1) - K20*A(2)

$ERROR
LOGADD = THETA(4)

IPRED = 0.000001 ;small number to prevent LOG(0)
IF(F.GT.0) IPRED = LOG(F)
W = LOGADD
IF(W.LE.1E-9) W=1E-9
Y = IPRED + W*EPS(1)
IRES = DV - IPRED
IWRES = IRES/W

$THETA
(0, 1) ;KA
(0, 1) ;VC
(0, 1) ;CL
(0, 0.1) ;LOGADDERR

$OMEGA
(0.1) ;IIV_KA
(0.1) ;IIV_VC
(0.1) ;IIV_CL

$SIGMA
1 FIX

Inter-Occasion Variability (IOV)

Inter-occasion variability reflects the variation of parameters between dosing occasions. Sometimes, the same subject, given the same dose of drug on two separate dosing occasions will have differing pharmacokinetic profiles. It may seem as though the amount of drug the subject was administered changed on each occasion! This may be encountered in studies where drugs are administered in very small volumes, or are subject to extremely nonlinear clearance mechanisms that magnify any effects of small variations in dosing. IOV, like IIV, is used to model these variations between occasions as being derived from a normal/log-normal distribution around a typical value.

Intuitively, IOV can present on parameters that may be expected to change dose-to-dose (perhaps due to slight variability in drug or drug administration). Parameters that frequently present with IOV are F_rel, k_a, D(N), R(N), etc. It is more difficult to justify, but by no means is it inappropriate, to include IOV terms on parameters generally associated with subject disposition, like Vc/F or CL/F. Some parameters, like CL/F, may vary from day-to-day due to liver status, diet, etc.

Note that IOV will require a column designating dosing occasions (OCC).

Inter-Occasion Variability on Relative Bioavailability

F_rel is a parameter reflecting the relative bioavailability and the typical value is FIXED to 1. The IOV on F_rel is log-normally distributed. We will assume there are three dosing occasions.

We use a clever trick in the \(\$\)OMEGA matrix to encode the IOV. Basically, the IOV ETAs should be sampled from the same distribution at each occasion. However, ETA sampling of \(\$\)OMEGAs only occurs once for each subject ID.

The solution is to assign three \(\$\)OMEGAs for the IOV variance, one for each OCC. One will be estimated, and the other two will be set to be the same as the estimated variance using the \(\$OMEGA~SAME\) option. Then, ON/OFF flags for each occasion are used to ensure only the occasion-specific value of ETA is being applied to F1. The flags are defined in the \(\$PK\) block.

$PROBLEM IOV on FREL

$INPUT IGNORE=C ID TIME AMT DV CMT EVID MDV OCC

$DATA nonmem_ready_dataset.csv

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(DEPOT)
COMP=(CENTRAL)
 
$PK
;OCC FLAGS
OCC1=0
OCC2=0
OCC3=0

IF(OCC.EQ.1)OCC1=1
IF(OCC.EQ.2)OCC2=1
IF(OCC.EQ.3)OCC3=1

;Typical Values
TVKA = THETA(1)
TVVC = THETA(2)
TVCL = THETA(3)
TVFREL = THETA(4)

;Parameters
KA = TVKA*EXP(ETA(1))
VC = TVVC*EXP(ETA(2))
CL = TVCL*EXP(ETA(3))

IOVF1 = OCC1*ETA(4) + OCC2*ETA(5) + OCC3*ETA(6)
F1 = TVFREL*EXP(IOVF1)

;Derived Parameters
K12 = KA
K20 = CL/VC

$DES
DADT(1) = -K12*A(1)
DADT(2) =  K12*A(1) - K20*A(2)

$ERROR
ADD = THETA(5)
PROP = THETA(6)

IPRED = F
W = PSQRT(ADD*ADD + PROP*PROP*IPRED*IPRED)
IF(W.LE.1E-9) W=1E-9
Y = IPRED + W*EPS(1)
IRES = DV - IPRED
IWRES = IRES/W

$THETA
(0, 1) ;KA
(0, 1) ;VC
(0, 1) ;CL
(1) FIXED ;FREL
(0, 0.1) ;ADDERR
(0, 0.1) ;PROPERR

$OMEGA
(0.1) ;IIV_KA
(0.1) ;IIV_VC
(0.1) ;IIV_CL

$OMEGA BLOCK(1) 0.1 ;IOV_F1_OCC1
$OMEGA BLOCK(1) SAME ;IOV_F1_OCC2
$OMEGA BLOCK(1) SAME ;IOV_F1_OCC3


$SIGMA
1 FIX

Posterior ETAs and IOV

Be cautious! If not every subject undergoes every dosing occasion in the dataset, the IOV ETA values associated with the missing occasions will collapse to 0 due to lack of information. This will influence the estimates of ETABAR and shrinkage. One can use the ETASTYPE=1 option in \(\$EST\) or the ETASXI(N) command to exclude specific IOV ETAs from subjects missing specific occasions.