Certara NLME Engine Knowledgebase

Knowledge Base Structure (for each entry):

Title: A concise title describing the topic.
Keywords: Relevant keywords for searching/retrieval.
Summary: A brief (1-2 sentence) summary of the topic.
Details: A more detailed explanation, including syntax, examples, and related concepts.
Related Topics: Links (using Titles) to other related entries in the knowledge base.
NONMEM Equivalent (Optional): If applicable, a brief description of the equivalent concept in NONMEM.

Part 1: Basic Model Structure and Statements

Title: Basic PML Model Structure

Keywords: model, structure, test, block, curly braces, main block

Summary: A PML model defines the structure of a pharmacokinetic/pharmacodynamic model. The basic structure consists of a main block enclosed in curly braces, typically named test(), containing statements that define the model’s components.

Details: The fundamental structure of a PML model is:

modelName() {
  // Statements defining the model
}

modelName(): The name of the model (e.g., test, oneCompartmentModel). The parentheses () are required. While any valid name can be used, test() is a common convention for simple models.
{ ... }: Curly braces define the model block. All model statements must be within this block.
Statements within the block define the model’s components (e.g., parameters, equations, error models). Statement order is generally not important, except within sequence blocks and for assignment statements.

Related Topics: Statements, Blocks, Comments, Variables

Title: Statements

Keywords: statement, assignment, declaration, semicolon, order

Summary: Statements are the building blocks of a PML model. They define parameters, equations, and other aspects of the model. Statements can span multiple lines, and semicolons are optional separators.

Details:

Types of Statements: PML includes various statement types, including:
- Assignment statements (e.g., C = A1 / V)
- Declaration statements (e.g., stparm, fixef, ranef, covariate, error)
- Control flow statements (within sequence blocks: if, while, sleep)
- Model component statements (e.g., deriv, observe, multi, LL, dosepoint, cfMicro)
Order: In general, the order of statements does not matter in PML, except:
- Within sequence blocks, where statements are executed sequentially.
- For assignment statements, where the order of assignments to the same variable matters.
- When defining micro constants before cfMicro
Semicolons: Semicolons (;) are optional separators between statements. They can improve readability but are not required.
Line Boundaries: Statements can span multiple lines.

Related Topics: Basic PML Model Structure, Blocks, Variables, Assignment Statements

Title: Blocks

Keywords: block, curly braces, sequence, scope

Summary: Blocks in PML are sets of statements enclosed in curly braces {}. They are used to group related statements and, in the case of sequence blocks, to define a specific execution order.

Details:

Main Model Block: The entire model is enclosed in a block (e.g., test() { ... }).
sequence Block: The sequence block is a special type of block that defines a sequence of statements to be executed in order, at specific points in time. This is used for handling discontinuous events and time-dependent actions.
Other Blocks: There are no other named blocks besides the main model block and sequence blocks.
dobefore, doafter blocks: used to define some actions related to dose delivery and observations.

Example (sequence block):

sequence {
  A = 10  // Initialize compartment A
  sleep(5) // Wait for 5 time units
  A = A + 5 // Add to compartment A
}

Related Topics: Basic PML Model Structure, Statements, sequence Block

Title: Comments

Keywords: comment, documentation, #, /* */, //

Summary: Comments are used to document the PML code, explaining its purpose and functionality. PML supports three comment styles: R-style, C-style, and C++-style.

Details:

R-style: # comment ... end-of-line (Everything after # on a line is a comment).
C-style: /* comment ... */ (Multi-line comments, cannot be nested).
C++-style: // comment ... end-of-line

Example:

test() {
  # This is an R-style comment
  /* This is a
     multi-line C-style comment */
  deriv(A = -k * A)  // C++-style comment
  fixef(k = c(,3,))
  x = 5 //valid statement
  y = 1 #valid statement
}

Important Note: While semicolons (;) are used as optional statement separators in PML, they are NOT reliable comment markers. Do not use a single semicolon to start a comment. The following examples demonstrate correct and incorrect comment usage:

test() {
       deriv(A1 = k)
       stparm(k = tvk) # This works because 'stparm' is valid.
       ; This is NOT a valid comment and will cause an error
       # This is a valid R-style comment.
       // This is a valid C++-style comment.
       \* This is valid C-style multi-line comment.
       */
       fixef(tvk = c(0,1,2))
}

Related Topics: Basic PML Model Structure

Title: Variables

Keywords: variable, declaration, assignment, scope, real, double, declared, functional

Summary: Variables in PML represent quantities within the model, such as parameters, compartment amounts, and concentrations. Variables can be declared explicitly or defined through assignment.

Details:

Data Types: All variables in PML are double-precision floating-point numbers (you can declare them using real or double, which are equivalent).
Declared Variables: Introduced by declaration statements like deriv, real, stparm, sequence, and covariate. These can be modified at specific points in time (e.g., within sequence blocks). Examples:
- deriv(A = -k * A) (Declares A as an integrator variable)
- real(x) (Declares x as a real variable, modifiable in sequence blocks)
- stparm(V = tvV * exp(nV)) (Declares V as a structural parameter)
Functional Variables: Introduced by assignment at the top level of the model (i.e., not within a sequence block). These are considered to be continuously calculated and cannot be modified within sequence blocks. Example:
- C = A / V (Defines C as the concentration, calculated from A and V)
Variable Names: Variable names must:
- Be case-sensitive
- Not contain special characters like “.”
- Not begin with an underscore “_”
Predefined variables PML contains predefined variables like t

Related Topics: Statements, Assignment Statements, Declaration Statements, stparm, fixef, ranef, deriv, real

Title: Predefined Variables

Keywords: predefined variables, t, time, built-in variables

Summary: PML has a limited set of predefined variables that are automatically available within the model.

Details:

t: Represents time. This variable is automatically available in time-based models (models containing at least one deriv statement).
The use of t on the right-hand-side of deriv statement is not recommended.
Important Note: Contrary to some other modeling languages, pi (the mathematical constant ≈ 3.14159) is NOT a predefined variable in PML. You must define it explicitly (e.g., pi = 3.141592653589793).

Related Topics: Variables, Time-Based Models, deriv Statement

Title: Assignment Statements

Keywords: assignment, =, equation, expression, continuous

Summary: Assignment statements are used to define the value of a variable based on an expression. They use the = operator. Variables defined by top-level assignment are considered continuously calculated.

Details:

Syntax: variable = expression
expression: Can be any valid mathematical expression involving numbers, variables, operators, and functions.
Top-Level Assignment: Assignments made at the top level of the model (outside of sequence blocks) define functional variables. These are continuously updated as the values of the variables they depend on change. Functional variables are analogous to variables defined by simple assignment within NONMEM’s $PK or $PRED blocks (before the Y= line).
Order matters: The order of the statements matter.
Multiple assignments: It is allowable to have multiple assignment statements assigned to the same variable, in which case the order between them matters.
Important Note: Variables defined and assigned via top-level assignment cannot be modified within sequence blocks.

Example:

test() {
  A = 10        # Assigns the value 10 to A
  V = 5         # Assigns the value 5 to V
  C = A / V     # Calculates C (concentration) continuously
  C = C + 2     # C will be equal to A/V + 2
}

Related Topics: Variables, Statements

Part 2: Parameter Declarations (stparm, fixef, ranef)

Title: stparm Statement

Keywords: structural parameter, stparm, group, parameter, population, fixed effect, random effect, covariate, IIV

Summary: The stparm statement defines structural parameters by linking the model’s structural level (compartments, equations) to the population level (fixed effects, random effects, and covariates).

Details:

A PML model consists of two levels: * The structural level, which contains compartments, doses, observations, and the core parameters describing the processes (structural parameters). * The population level, which includes typical values of parameters (fixed effects), inter-individual variation (random effects), and covariate effects.

The stparm statement is the statement that ties these two levels together.

Purpose: Defines structural parameters and their relationships to fixed effects, random effects, and covariates for estimation and simulation.
Syntax: stparm(parameter = expression)
- parameter: The name of the structural parameter (e.g., V, Cl, Ka, F).
- expression: An expression defining the parameter, which can include:
  - Fixed effects (typically tv prefix, e.g., tvV).
  - Random effects (typically n prefix, e.g., nV).
  - Covariates and group parameters.
  - Mathematical operators and functions.
Execution and Performance: The stparm statements are only re-executed during an iteration when a covariate value changes. They are not evaluated at every step of the ODE solver. Therefore, they should not include variables that need to be evaluated continuously, like concentration (C).

The `group` Statement

To simplify stparm statements and improve performance, you can pre-calculate parts of them in a group block.

Purpose: A group statement defines group parameters that can be used in stparm expressions. It ensures these calculations are done efficiently.
Execution: group statements are executed before stparm statements and are also only re-evaluated when a covariate changes.
Constraints: Expressions within group can only be functions of covariates and fixed effects.
Example (Performance): Calculating Body Surface Area (BSA) should be done in a group statement to avoid repeated calculations within the ODE solver.

group(BSA = (WT^0.5378) * (HT^0.3964) * 0.024265)     
stparm(Cl = tvCl * (BSA/1.73)^0.75 * exp(nCl))

Example (Clarity):

group {       
  tvVMale   = tvV
  tvVFemale = tvV * dVdGender     
}     

stparm(V = ((Gender == 0) ? tvVMale : tvVFemale) * exp(nV))

Common Parameterization Styles

Log-Normal: parameter = tvParameter * exp(nParameter) (For positive parameters like Cl, V).
Logit-Normal: parameter = ilogit(tvParameter + nParameter) (Recommended for parameters constrained between 0 and 1, like bioavailability F).
Normal: parameter = tvParameter + nParameter (For parameters that can be positive or negative, like E0).
Multiple stparm Statements: A model can have multiple stparm statements.
Important Note (NONMEM): stparm defines estimable structural parameters. A simple variable assignment in NONMEM’s $PK block that is not a THETA should be a top-level assignment in PML (e.g., C = A1/V), not a stparm statement.
Conditional Logic: Use the ternary operator (? :) for conditional logic within stparm expressions, not if/else.
Important Note (Time-Varying): When a stparm depends on time t or an interpolate() covariate, closed-form solutions are not applicable, and the model will rely on a numerical ODE solver.

Example:

stparm(V  = tvV * exp(nV))                  // Volume (log-normal) 
stparm(Ka = tvKa)                            // Absorption rate constant (fixed effect) 
stparm(Cl = tvCl * exp(dCldSex1*(Sex==1)) * exp(nCl + nClx0*(Period==1) + nClx1*(Period==2))) // Complex clearance

Related Topics: group Statement, fixef Statement, ranef Statement, Covariates, Fixed Effects, Random Effects, Log-Normal Distribution, if and Ternary Operator

Title: `fixef` Statement

Keywords: fixed effect, fixef, population parameter, typical value, initial estimate, bounds, enable

Summary: The fixef statement defines fixed-effect parameters (often called “thetas” in NONMEM), representing the typical or population values of structural parameters. It also allows specifying initial estimates, bounds, and enabling/disabling fixed effects for covariate search.

Details:

Purpose: Defines fixed-effect parameters and their properties.
Syntax: fixef(parameter[(freeze)][(enable=c(int))] [= c([lower bound],[ initial estimate],[ upper bound])])
- parameter: The name of the fixed-effect parameter (e.g., tvV, tvCl, dCldSex).
- freeze: (Optional) If present, the parameter is not estimated; it’s fixed at its initial estimate. May use freeze for parameters that are defined and assigned values in the NONMEM $PK block but are not THETAs (i.e., not estimated), but probably it is better to use the direct assignment outside fixef/stparm statements
- enable=c(int): (Optional) Controls whether the parameter is considered during covariate search procedures (0 = disabled, 1 = enabled). Has no effect on a regular model fit.
- = c([lower bound],[ initial estimate],[ upper bound])]: (Optional) Specifies the lower bound, initial estimate, and upper bound. Any of these can be omitted (using commas as placeholders). If only one value is given without bounds, it represents the initial value
Default Initial Values: If no initial estimate is provided:
- Covariate effects: Default initial value is 0.
- Other fixed effects: Default initial value is 1.
Multiple fixef: It is possible to have more than one fixef statement.
Critical Restriction: Use Numeric Literals Only The values provided for the lower bound, initial estimate, and upper bound within a fixef statement must be numeric literals (i.e., plain numbers). You cannot use function calls (like log(), sqrt()), variables, or mathematical expressions. If you are re-parameterizing a model (e.g., estimating a log-transformed parameter), you must perform the calculation before writing the PML code and insert the resulting number.
Example:

# GOAL: Re-parameterize tvMTT and provide an initial estimate of log(3.27)     
# INCORRECT - Function call inside fixef() will cause an error     
stparm(MTT = exp(tvLogMTT))     
fixef(tvLogMTT = c(, log(3.27), ))     
# CORRECT - Pre-calculate the value and use the number directly     
# Pre-calculation: log(3.27) is approximately 1.1848     
stparm(MTT = exp(tvLogMTT))     
fixef(tvLogMTT = c(, 1.1848, ))

Example:

    fixef(tvV  = c(0, 10, 100))  // Initial estimate 10, bounds [0, 100]
    fixef(tvCl = c(, 5, )) // Initial estimate 5, no bounds
    fixef(dCldSex1(enable=c(0)) = c(, 0, )) // Enabled for covariate search, initial estimate 0

Related Topics: stparm Statement, ranef Statement, Covariates, Fixed Effects, Bounds, Covariate Search

Title: `ranef` Statement

Keywords: random effect, ranef, inter-individual variability, IIV, variance, covariance, omega, diag, block, same

Summary: The ranef statement defines random-effect parameters (often called “etas” in NONMEM), representing inter-individual variability (IIV) in structural parameters. It also defines the variance-covariance structure of the random effects.

Details:

Purpose: Defines random effects and their variance-covariance matrix.
Syntax: ranef(specification)
- specification: Defines the random effects and their relationships. Common options include:
  - diag(parameter1, parameter2, ...): Diagonal covariance matrix (random effects are independent).
  - block(parameter1, parameter2, ...): Full covariance matrix (random effects can be correlated).
  - same(parameter1, parameter2, ...): Specifies that parameters share the same variance (and covariance, if within a block).
  - = c(...): initial estimates for variance/covariance
The initial estimates are provided for variance-covariance matrix
Multiple ranef Statements: A model can have multiple ranef statements.
NONMEM Equivalent:
- diag(nV, nCl) is similar to NONMEM’s:
```
$OMEGA DIAG
0.04 ;nV
0.09 ;nCl
```
- block(nV, nCl) is similar to NONMEM’s:
```
$OMEGA BLOCK(2)
0.04 ;nV
0.02 0.09 ;nCl
```

Examples:

ranef(diag(nV, nCl) = c(0.04, 0.09))  // Diagonal covariance: V and Cl vary independently

ranef(block(nV, nCl) = c(0.04, 0.02, 0.09))  // Full covariance: V and Cl are correlated

ranef(diag(nV) = c(0.04), diag(nCl) = c(0.09)) // Equivalent to the first example

ranef(block(nCl1, nCl2) = c(1, 0.5, 2), same(nCl3, nCl4)) //block + same

Related Topics: stparm Statement, fixef Statement, Random Effects, Variance, Covariance, Inter-Individual Variability

Title: Fixed Effects

Keywords: fixed effect, population parameter, typical value, theta, fixef

Summary: Fixed effects represent the typical or population values of parameters in the model. They are estimated from the data and are assumed to be constant across all individuals.

Details: * Represented by fixef * Usually defined with tv prefix

Related Topics: fixef Statement, stparm Statement, Random Effects

Title: Random Effects

Keywords: random effect, inter-individual variability, IIV, eta, ranef, variance, covariance

Summary: Random effects represent the variability in parameters between individuals. They are assumed to be drawn from a distribution (usually normal or log-normal) with a mean of zero and an estimated variance-covariance matrix.

Details: * Represented by ranef * Usually defined with n prefix

Related Topics: ranef Statement, stparm Statement, Fixed Effects, Variance, Covariance, Inter-Individual Variability

Title: Log-Normal Distribution (in Parameter Definitions)

Keywords: log-normal, distribution, positive parameter, stparm, exp

Summary: The log-normal distribution is commonly used for parameters that must be positive (e.g., clearance, volume). It’s implemented in stparm using the exp() function.

Details:

Syntax: parameter = tvParameter * exp(nParameter)
- tvParameter: The typical value (fixed effect).
- nParameter: The random effect (normally distributed with mean 0).
- exp(): The exponential function. This ensures that parameter is always positive, regardless of the value of nParameter.
Why Log-Normal? Many pharmacokinetic parameters (e.g., clearance, volume) can only have positive values. The log-normal distribution guarantees positivity.

Example:

stparm(
  Cl = tvCl * exp(nCl)  // Clearance is log-normally distributed
)

Related Topics: stparm Statement, Fixed Effects, Random Effects, Normal Distribution

Title: Normal Distribution (in Parameter Definitions)

Keywords: normal, distribution, stparm, additive

Summary: The normal distribution can be used for parameters that can take on both positive and negative values.

Details:

Syntax: parameter = tvParameter + nParameter

Example:

stparm(
  Effect = tvEffect + nEffect  // Effect can be positive or negative
)

Related Topics: stparm statement, Log-Normal Distribution

Title: Bounds (in `fixef` statements)

Keywords: bounds, fixef, lower bound, upper bound, parameter constraints

Summary: Bounds (lower and upper limits) can be specified for fixed effects within the fixef statement. This constrains the parameter values during estimation, preventing them from taking on unrealistic values.

Details:

Syntax: fixef(parameter = c([lower bound],[ initial estimate],[ upper bound]))
Any value could be skipped

Example:

fixef(
  tvCl = c(0, 5, 20)  // Cl must be between 0 and 20, with an initial estimate of 5
)

Related Topics: fixef statement

Part 3: Covariates

Title: Covariates

Keywords: covariate, independent variable, predictor, time-varying, categorical, continuous, occasional, fcovariate, covariate, input, data mapping

Summary: Covariates are variables that can influence the model’s parameters or outcomes. They represent characteristics of the individuals, their environment, or external influences. PML supports continuous, categorical, and occasional covariates, any of which can be time-varying.

Details:

Purpose: To incorporate the effects of individual characteristics, external factors, or time-dependent changes on the model.
Types of Covariates:
- Continuous: Take on a continuous range of numerical values (e.g., weight, age, creatinine clearance).
- Categorical: Take on a limited number of discrete values representing categories (e.g., sex, disease status, treatment group).
- Occasional: Represent different occasions or periods within an individual’s data (e.g., different treatment cycles). These are inherently time-varying.
Declaration: Covariates must be declared using either the covariate, fcovariate, or interpolate statement.
Time-Varying Behavior: Any type of covariate (continuous, categorical, or occasional) can be time-varying. The covariate, fcovariate, and interpolate statements control how the covariate values are handled over time.
Usage: Covariates can be used in:
- stparm statements to model their effects on structural parameters.
- Expressions within the model (e.g., in deriv statements or likelihood calculations).
Data Mapping: Covariates are linked to columns in the input dataset through column mappings. This is typically done in a separate file or within the user interface of the modeling software (e.g., Phoenix NLME). The correct syntax for mapping is, for example, covr(wt<-"wt").

Related Topics: covariate Statement, fcovariate Statement, interpolate Statement, stparm Statement, Continuous Covariates, Categorical Covariates, Occasional Covariates, Time-Varying Covariates, Data Mapping

Title: `covariate` Statement

Keywords: covariate, backward extrapolation, time-varying covariate

Summary: The covariate statement declares a covariate and specifies that its values should be extrapolated backward in time. The most recent value of the covariate is used until a new value is encountered.

Details:

Syntax: covariate(covariateName) or covariate(covariateName())
- covariateName: The name of the covariate.
- (): Empty parentheses are required for categorical covariates and recommended for occasional covariates (though technically optional for occasion).
Backward Extrapolation: The value of the covariate at any given time is the most recent value observed before that time. If no value is given at time=0, the first available value is used.

Example:

covariate(Weight)  // Declares 'Weight' as a continuous covariate
covariate(Sex())    // Declares 'Sex' as a categorical covariate

NONMEM Equivalent: There is no direct equivalent in NONMEM. NONMEM’s default behavior for covariates is backward extrapolation, similar to PML’s covariate.

Related Topics: Covariates, fcovariate Statement, interpolate Statement, Time-Varying Covariates, Categorical Covariates, Occasional Covariates

Title: `fcovariate` Statement

Keywords: fcovariate, forward extrapolation, time-varying covariate, occasion covariate

Summary: The fcovariate statement declares a covariate and specifies that its values should be extrapolated forward in time. The current value is used until a new value is encountered. fcovariate is generally preferred for occasion covariates.

Details:

Syntax: fcovariate(covariateName) or fcovariate(covariateName())
- covariateName: The name of the covariate.
- (): Empty parentheses are required for categorical covariates and recommended (but technically optional) for occasional covariates.
Forward Extrapolation: The value of the covariate at any given time is the value observed at that time, and this value is carried forward until a new value is encountered. The first available value is also carried backward if no value is given at time=0.
Occasion Covariates: fcovariate is generally preferred for occasion covariates.

Example:

fcovariate(DoseRate)   // Declares 'DoseRate' as a time-varying covariate
fcovariate(Occasion())  // Declares 'Occasion' as an occasion covariate
fcovariate(Occasion) // also valid, but less explicit

NONMEM Equivalent: There’s no direct equivalent in NONMEM. You would typically handle forward extrapolation of occasion covariates implicitly through the structure of your dataset and control stream.

Related Topics: Covariates, covariate Statement, interpolate Statement, Time-Varying Covariates, Occasional Covariates, Categorical Covariates

Title: `interpolate` Statement

Keywords: interpolate, covariate, linear interpolation, time-varying covariate, continuous covariate

Summary: The interpolate statement declares a continuous covariate whose values are linearly interpolated between the time points at which it is set.

Details:

Syntax: interpolate(covariateName)
- covariateName: The name of the covariate.
Linear Interpolation: The value of the covariate varies linearly between the time points at which it is set in time-based models.
Extrapolation: If no covariate value is given, the latest value is carried forward. If no value is given at time=0, the first available value is used.
Continuous Covariates Only: interpolate can only be used with continuous covariates.

Example:

interpolate(InfusionRate)

NONMEM Equivalent: There is no direct equivalent in NONMEM. Linear interpolation of covariates is not a built-in feature. You would typically pre-process your data to create interpolated values if needed.

Related Topics: Covariates, covariate Statement, fcovariate Statement, Time-Varying Covariates, Continuous Covariates

Title: Continuous Covariates

Keywords: continuous covariate, covariate, numerical, time-varying

Summary: Continuous covariates take on a continuous range of numerical values. They can be time-varying or constant within an individual.

Details:

Examples: Weight, age, creatinine clearance, body surface area.
Declaration: Declared using covariate, fcovariate, or interpolate.
Usage: Used directly in mathematical expressions within stparm statements or other model equations.
Mapping: covr(Wt<-"Wt")

Example:

fcovariate(Weight)  // Weight as a time-varying covariate

stparm(
  Cl = tvCl * (Weight / 70)^0.75 * exp(nCl)  // Allometric scaling of clearance
)

Related Topics: Covariates, covariate Statement, fcovariate Statement, interpolate Statement, stparm Statement, Data Mapping

Title: Categorical Covariates

Keywords: categorical covariate, covariate, discrete, categories, dummy variable, indicator variable, time-varying, data mapping

Summary: Categorical covariates take on a limited number of discrete values representing categories. PML requires you to define the mapping between the data values and category names in the column definition file.

Details:

Examples: Sex (Male/Female), Disease Status (Normal/Mild/Severe), Treatment Group (Placebo/Active).
Declaration: Declared using covariate(CovariateName()) or fcovariate(CovariateName()). The empty parentheses () are required.
Mapping (Column Definition File): You must map the covariate and its categories in the column definition file (or equivalent interface in your software).
- With Labels: If your data file already contains meaningful category labels (e.g., “Male”, “Female”), map them directly. The general NONMEM-style syntax is: covr(Sex <- "Sex"("Male" = 0, "Female" = 1)) // Example This maps the “Sex” column in the data to a covariate named Sex in the model, with “Male” coded as 0 and “Female” as 1. The first category is used as a reference category.
- Without Labels: If your data file uses numeric codes (e.g., 0, 1) without explicit labels, you can define the labels during mapping using empty parentheses: covr(Sex <- "Sex"()) // Example In that case the first unique value will be used as a reference.
Usage in stparm: You typically use logical expressions (CovariateName == value) within stparm statements to model the effects of different categories. This creates implicit “dummy variables.” The first category encountered in the data is treated as the reference category, and fixed effects for other categories represent differences from the reference. Use the ternary operator (? :) for more complex conditional logic.

Example (PML code):

test() {
  fcovariate(Sex())  // Declare Sex as a categorical covariate

  stparm(
    Cl = tvCl * exp(dClSex * (Sex == 1) + nCl)  // Effect of Sex on Cl
  )

  fixef(
    tvCl   = c(, 10, ),
    dClSex = c(, 0, )  // Fixed effect for Sex=1 (relative to Sex=0)
  )

  ranef(diag(nCl) = c(0.25))

  // ... rest of the model ...
}

NONMEM Equivalent: The PML code above is similar to the following in NONMEM (using abbreviated code):

$INPUT ... SEX ...
$SUBROUTINES ...
$PK
  TVCL = THETA(1)
  DCLSEX = THETA(2)
  CL = TVCL * EXP(DCLSEX*(SEX.EQ.1) + ETA(1))
$THETA
 (0, 10) ; TVCL
 (0)      ; DCLSEX
$OMEGA
 0.25 ; ETA(1) variance (IIV on CL)

Related Topics: Covariates, covariate Statement, fcovariate Statement, stparm Statement, Fixed Effects, Data Mapping, Dummy Variables, if and Ternary Operator

Title: Occasion Covariates and Inter-Occasion Variability (IOV)

Keywords: occasion covariate, IOV, inter-occasion variability, fcovariate, random effects, fixed effects, within-subject variability

Summary: Occasion covariates represent distinct periods in a subject’s data (e.g., treatment cycles). Their primary use is to model Inter-Occasion Variability (IOV), which is the random, unpredictable variability within a single subject from one occasion to the next. This is modeled using random effects, not fixed effects.

Details:

Concept: While Inter-Individual Variability (IIV) describes why Subject A’s average parameter value is different from Subject B’s, IOV describes why Subject A’s parameter value during Occasion 1 is different from their own value during Occasion 2.
Declaration: Always declare the occasion variable, preferably with fcovariate(Occasion()). The parentheses () are recommended for clarity.

Correct Implementation: Modeling IOV with Random Effects

This is the standard and pharmacokinetically correct approach. It assumes that a subject’s parameter value for a specific occasion is a random deviation from their overall mean parameter value.

PML Syntax: You create a separate random effect (eta) for each occasion level and add it to the structural parameter expression.
- stparm(Param = tvParam * exp( nParam + nParamOcc1*(Occasion==1) + nParamOcc2*(Occasion==2) + ... ))
The term nParam represents the subject’s IIV (their deviation from the population typical value tvParam).
The terms nParamOcc1, nParamOcc2, etc., represent the IOV (the deviation for that specific occasion from the subject’s mean).

Example (IOV on Clearance):

test() {
  fcovariate(Occasion()) // Assume Occasion has levels 1, 2, 3

  // Cl includes a random effect for IIV (nCl) and separate random effects for IOV on each occasion
  stparm(Cl = tvCl * exp( nCl + nCl_Occ1*(Occasion==1) + nCl_Occ2*(Occasion==2) + nCl_Occ3*(Occasion==3) ))

  fixef(tvCl = c(, 1, ))
  ranef(diag(nCl) = c(1)) // IIV variance on Cl

  // Define the IOV variances. 'same()' is often used to assume variability is equal across occasions.
  ranef(diag(nCl_Occ1) = c(1), same(nCl_Occ2), same(nCl_Occ3))
  ...
}

Incorrect Implementation: Modeling Occasions with Fixed Effects

It is also possible to model an occasion as a fixed effect, but this answers a different scientific question and is generally not what is meant by an “occasion effect.” A fixed effect tests whether the population average parameter is systematically different on one occasion versus another (e.g., “Is clearance for all subjects 20% lower on Occasion 2?”). This is a cohort effect, not within-subject variability.

Example of an (often incorrect) fixed-effect model:

// This model tests if the POPULATION mean Cl 
// is different on Occasion 2 vs Occasion 1
stparm(Cl = tvCl * exp( dCldOcc2*(Occasion==2) ) * exp(nCl))
fixef(dCldOcc2 = c(, 0, )) // Fixed effect for Occasion 2

Related Topics: fcovariate Statement, Random Effects, ranef Statement, Fixed Effects, stparm Statement, Inter-Individual Variability

Title: Time-Varying Covariates

Keywords: time-varying covariate, covariate, fcovariate, covariate, changing covariate, dynamic covariate

Summary: Time-varying covariates are covariates whose values can change within an individual over time. They are essential for modeling dynamic processes and external influences that vary during the observation period. Any type of covariate (continuous, categorical, or occasional) can be time-varying.

Details:

Declaration: Declared using covariate, fcovariate, or interpolate. fcovariate is generally preferred for most time-varying covariates, especially for categorical and occasion covariates.
Data: The input data must include multiple records per individual, with different time points and corresponding covariate values.
Extrapolation:
- covariate: Backward extrapolation.
- fcovariate: Forward extrapolation (and backward extrapolation for the first value).
- interpolate: Linear interpolation between defined points (continuous covariates only), forward extrapolation after the last defined point.

Example:

fcovariate(DoseRate)   // DoseRate can change over time

stparm(
  Cl = tvCl * exp(dClDoseRate * DoseRate + nCl)  // Clearance depends on DoseRate
)

// ... rest of the model ...

Related Topics: Covariates, covariate Statement, fcovariate Statement, interpolate Statement, Continuous Covariates, Categorical Covariates, Occasional Covariates

Part 4: Structural Model Definition

Title: `deriv` Statement

Keywords: deriv, differential equation, ODE, integrator, state variable, time-based model, dynamic model

Summary: The deriv statement defines an ordinary differential equation (ODE) in the model. This is used to describe how the amount of drug (or other quantities) in a compartment changes over time. Models with deriv statements are inherently time-based.

Details:

Purpose: To define the rate of change of a state variable (typically a compartment amount).
Syntax: deriv(variable = expression)
- variable: The name of the state variable being differentiated (e.g., A1, Aa, CumHaz). This variable is often called an “integrator variable.”
- expression: An expression defining the rate of change of the variable. This can involve parameters, covariates, other variables, and mathematical functions.
Time-Based Models: If a model contains any deriv statement, it is considered a time-based model, and the built-in variable t (representing time) is automatically available.
State variable modification Variables on the left side of deriv statements can be modified when the model is stopped
Multiple deriv statements: A model will typically have multiple deriv statements, one for each compartment or state variable being modeled.
Right-hand-side Discontinuity: The use of t variable on the right-hand-side is discouraged.

Example:

deriv(Aa = -Ka * Aa)      // Rate of change of amount in absorption compartment Aa
deriv(A1 = Ka*Aa -Cl * C)       // Rate of change of amount in compartment A1

NONMEM Equivalent: The PML code above is similar to the following NONMEM code:

$DES
DADT(1) = -KA*A(1)
DADT(2) = KA*A(1) -CL*A(2)/V

Related Topics: Time-Based Models, Compartment Models, dosepoint Statement, urinecpt Statement, State Variables

Title: Compartment Models

Keywords: compartment model, compartment, differential equation, PK model, ADVAN

Summary: Compartment models are a common type of pharmacokinetic model where the body is represented as a series of interconnected compartments. Drug movement between compartments is described by differential equations.

Details:

Compartments: Represent theoretical spaces where the drug distributes (e.g., central compartment, peripheral compartment, absorption compartment).
Differential Equations: Describe the rate of change of the amount of drug in each compartment.
Parameters: Define the rates of transfer between compartments and elimination from the body (e.g., clearance, volume, rate constants).
Implementation in PML: Compartment models can be implemented using:
- deriv statements (for custom models or when flexibility is needed).
- cfMicro or cfMacro statements (for standard 1, 2, or 3-compartment models).

Related Topics: deriv Statement, cfMicro Statement, cfMacro Statement, Time-Based Models, Pharmacokinetics

Title: `cfMicro` Statement

Keywords: cfMicro, closed-form solution, compartment model, micro-constants, deriv, time-varying covariate, ODE solver

Summary: Provides an efficient closed-form (analytical) solution for standard compartment models using micro-constants. This statement is not compatible with time-varying parameters and must be replaced by deriv statements in such cases.

Details:

Purpose: To define a standard one-, two-, or three-compartment model without explicitly writing out the differential equations, which can improve performance.
Parameterization: Uses micro-constants (e.g., Ke, K12, K21).
Syntax: cfMicro(A, Ke, [K12, K21], [K13, K31], [first = (Aa = Ka)])
- A: The amount in the central compartment.
- Ke: The elimination rate constant.
- K12, K21: (Optional) Transfer rate constants for a two-compartment model.
- K13, K31: (Optional) Transfer rate constants for a three-compartment model.
- first = (Aa = Ka): (Optional) Specifies first-order absorption from a depot compartment Aa with absorption rate constant Ka.

CRITICAL LIMITATION: Time-Varying Parameters

The cfMicro statement uses a piecewise closed-form (analytical) solution. The fundamental assumption is that model parameters (like Cl, V, K12, etc.) are constant between time-based events (such as doses or changes in covariate values). When an event occurs, the system re-evaluates the parameters and applies the closed-form solution for the next interval.

This means cfMicro can correctly handle structural parameters that depend on time-varying covariates declared with covariate or fcovariate, as these values change at discrete time points.

For example, a model with stparm(Cl = tvCl * (CRCL/100) * exp(nCl)) and covariate(CRCL) will work correctly with the following data, as Cl is constant between t=0, t=2, and t=4:

ID	TIME	AMT	Cobs	CRCL
1	0	100	0	70
1	2	.	80	70
1	4	.	70	60

However, this piecewise assumption is violated if a structural parameter changes continuously over time. This happens in two main scenarios:

A structural parameter depends on a covariate declared with interpolate().
A structural parameter’s definition explicitly includes the time variable t.

CONSEQUENCE: Using cfMicro with continuously time-varying parameters will produce incorrect results without generating a syntax error.
REQUIRED ACTION: If any model parameter changes continuously with time (due to interpolate() or the use of t), you MUST use deriv statements to define the model structure. This forces the use of a numerical ODE solver that can correctly handle the dynamic parameters.
Stricter Rules for cfMacro and cfMacro1: These are pure closed-form solutions and are more restrictive. They require structural parameters to be constant for the entire observation interval and do not support time-varying covariates of any kind.

Example of an INVALID model (due to interpolate):

test() {
    interpolate(scr) // Makes scr a CONTINUOUSLY time-varying covariate
    stparm(Cl = tvCl + dCldscr * scr) // Cl is now continuously time-varying
    cfMicro(A1, Cl / V) // INVALID: cfMicro cannot handle a continuously time-varying Cl
    ...
}

Example of the CORRECT equivalent model:

test() {
    interpolate(scr) // Continuously time-varying covariate
    stparm(Cl = tvCl + dCldscr * scr) // Cl is continuously time-varying
    deriv(A1 = - (Cl/V) * A1) // CORRECT: Must use deriv statement
    ...
}

pyDarwin Automation Note: When designing a model search, if any candidate model in your search space includes a time-varying interpolated covariate effect on a structural parameter, all structural model options must be written using deriv statements. Do not mix cfMicro and deriv based structural models in the same search if time-varying covariates are involved, as it would be an invalid comparison.

Related Topics: deriv Statement, Compartment Models, covariate Statement, interpolate Statement, fcovariate Statement, Time-Varying Covariates

Title: `cfMacro` and `cfMacro1` Statements

Keywords: cfMacro, cfMacro1, closed-form solution, compartment model, macro-constants, A, Alpha, B, Beta, C, Gamma, stripping dose

Summary: The cfMacro and cfMacro1 statements provide closed-form solutions for standard compartment models, parameterized using macro-constants (coefficients and exponents of exponential terms). cfMacro offers more options, including a stripping dose.

Details:

Purpose: To define a standard compartment model without explicitly writing out the differential equations, using a macro-constant parameterization.
cfMacro Syntax: cfMacro(A, C, DoseVar, A, Alpha, [B, Beta], [C, Gamma], [strip=StripVar], [first=(Aa=Ka)])
- A: The amount in the central compartment (cannot be referred to elsewhere in the model).
- C: The concentration variable in the model.
- DoseVar: A variable to record the initial bolus dose (needed for idosevar).
- A, Alpha: Macro-constants for the first exponential term.
- B, Beta: (Optional) Macro-constants for the second exponential term (two-compartment model).
- C, Gamma: (Optional) Macro-constants for the third exponential term (three-compartment model).
- strip=StripVar: (Optional) Specifies a covariate (StripVar) to provide a “stripping dose” for simulations.
- first = (Aa = Ka): (Optional) Specifies the first order absorption
cfMacro1 Syntax: cfMacro1(A, Alpha, [B, Beta], [C, Gamma], [first=(Aa=Ka)])
- Simplified version of cfMacro.
- A: Amount in the central compartment.
- Alpha, B, Beta, C, Gamma: Macro-constants. The response to bolus dose is predefined
- first = (Aa = Ka): (Optional) Specifies the first order absorption
Stripping Dose: The strip option in cfMacro allows you to specify a covariate that provides a “stripping dose” value. This is used in simulations to represent the initial dose used when the model was originally fit.

Example (cfMacro with two compartments and stripping dose):

cfMacro(A1, C1, A1Dose, A, Alpha, B, Beta, strip=A1Strip)
dosepoint(A1, idosevar = A1Dose)
covariate(A1Strip)

Example (cfMacro1 with one compartment):

cfMacro1(A, Alpha)

Related Topics: Compartment Models, deriv Statement, cfMicro Statement, Macro-Constants, Closed-Form Solutions, Pharmacokinetics

Title: Micro-Constants vs. Macro-Constants

Keywords: micro-constants, macro-constants, Ke, K12, K21, A, Alpha, B, Beta, C, Gamma, parameterization, cfMicro, cfMacro

Summary: Compartment models can be parameterized using either micro-constants (rate constants) or macro-constants (coefficients and exponents of exponential terms). cfMicro uses micro-constants, while cfMacro and cfMacro1 use macro-constants.

Details:

Micro-Constants:
- Rate constants that describe the transfer of drug between compartments and elimination from the body (e.g., Ke, K12, K21).
- More directly related to the underlying physiological processes.
- Used with cfMicro and deriv statements.
Macro-Constants:
- Coefficients and exponents of exponential terms in the closed-form solution (e.g., A, Alpha, B, Beta, C, Gamma).
- Less intuitive from a physiological perspective.
- Used with cfMacro and cfMacro1.
Conversion: It’s possible to convert between micro-constants and macro-constants, but the equations can be complex, especially for models with more than one compartment.

Related Topics: cfMicro Statement, cfMacro Statement, cfMacro1 Statement, Compartment Models, Pharmacokinetics

Title: dosepoint Statement

Keywords: dosepoint, dose, dosing, bolus, infusion, tlag, duration, rate, bioavail, dobefore, doafter, idosevar, infdosevar, infratevar, split, dosepoint1, dosepoint2

Summary: The dosepoint statement specifies a compartment that can receive doses (either bolus or infusion). It also allows for options like lag time, infusion duration/rate, bioavailability, and actions performed before/after dosing. Multiple dosepoint statements can be used to model complex administration schemes.

Details:

Purpose: To define where and how doses are administered in the model.
Syntax: dosepoint(compartmentName [, tlag = expr][, duration = expr][, rate = expr][, bioavail = expr][, dobefore = sequenceStmt][, doafter = sequenceStmt][, split][, idosevar = variableName][, infdosevar = variableName][, infratevar = variableName])
- compartmentName: The name of the compartment receiving the dose.
- tlag = expr: (Optional) Time delay before dose delivery.
- duration = expr: (Optional) Duration of a zero-order infusion.
- rate = expr: (Optional) Rate of a zero-order infusion. Cannot be used with duration.
- bioavail = expr: (Optional) Bioavailability fraction (0 to 1).
- idosevar = variableName: (Optional) Captures the value of the first bolus dose.
- infdosevar = variableName: Captures the value of the first infusion dose.
- infratevar = variableName: (Optional) Captures the infusion rate of the first infusion dose.
Bolus vs. Infusion:
- If neither duration nor rate is specified, the dose is treated as a bolus (instantaneous).
- If duration or rate is specified, the dose is treated as a zero-order infusion.

Modeling Multiple Administration Routes

PML provides powerful ways to model scenarios with more than one route of administration using dosepoint, dosepoint1 (a direct synonym), and dosepoint2. The correct approach depends on whether you are modeling mutually exclusive formulations or partitioning a single dose.

1. Modeling Independent Dose Streams (Mutually Exclusive Formulations)

This is the standard approach for modeling different, independent routes of administration, such as IV vs. Tablet vs. Capsule.

Concept: Use multiple dosepoint or dosepoint2 statements. Each statement is mapped to a separate dose column in the input dataset. This creates independent “streams” for each formulation.
dosepoint vs. dosepoint2: You can have multiple dosepoint statements as long as they target different compartments. If you need to dose into the same compartment from different streams (e.g., a Tablet and a Capsule both absorbing into Aa), you use dosepoint for the first stream and dosepoint2 for the second.

Example (IV, Tablet, Capsule):

    dosepoint(A1)                   // IV stream
    dosepoint(Aa, bioavail=F_tab)   // Tablet stream
    dosepoint2(Aa, bioavail=F_cap)  // Capsule stream
    ```
**Corresponding Column Mapping:**

dose(A1<-"DOSEIV")
dose(Aa<-"DOSETABLET")
dose2(Aa<-"DOSECAPSULE")
```

2. Splitting a Single Dose Amount (`split` argument)

This pattern is used for a single formulation that has multiple simultaneous absorption pathways (e.g., a pill with both immediate-release and extended-release components).

Concept: This is achieved by using dosepoint and dosepoint2 on the same compartment and adding the crucial split argument to the first dosepoint statement.
split Argument: This tells the system that the dose() and dose2() mappings in the column definition file will both point to the same data column (e.g., AMT). The bioavail option in each statement then determines what fraction of that single AMT value goes to each pathway.

Example (50% Bolus, 50% Infusion from one dose):

    // 50% of AMT is a bolus into Aa. 'split' is essential.
    dosepoint(Aa, bioavail=0.5, split)
    // The other 50% of AMT is an infusion into Aa.
    dosepoint2(Aa, bioavail=0.5, duration=3)
    ```
**Corresponding Column Mapping:**

dose(Aa<-"AMT", split) // Both map to the same "AMT" column
dose2(Aa<-"AMT")
```

Important Note: A dosepoint statement is REQUIRED for any compartment that receives doses. Without it, AMT values in the dataset for that compartment will be ignored.

NONMEM Equivalent: Handling multiple dose streams in NONMEM typically involves using the CMT data item and IF/ELSE logic in the control stream to direct doses, or by structuring the dataset with different EVID values. The split functionality has no direct equivalent and requires more complex data manipulation in NONMEM.

Related Topics: Dosing, Bolus Dose, Infusion, Bioavailability, sequence Block, Compartment Models, Modeling Complex Absorption Schemes

Title: `urinecpt` Statement

Keywords: urinecpt, elimination compartment, excretion, steady-state

Summary: The urinecpt statement defines an elimination compartment, similar to deriv, but with specific behavior during steady-state dosing: it’s ignored.

Details:

Purpose: To model the elimination of drug or a metabolite into an accumulating compartment (typically urine).
Syntax: urinecpt(variable = expression [, fe = fractionExp])
- variable: The name of the elimination compartment amount.
- expression: Defines the rate of elimination into the compartment.
- fe = fractionExp: (Optional) Specifies fraction excreted
Steady-State Behavior: The key difference between urinecpt and deriv is that during steady-state simulations, the urinecpt statement is ignored. This is because, at steady state, the amount in the elimination compartment is not relevant to the dynamics of the system.
Resetting: The amount in the compartment could be set to 0 after being observed using observe statement and doafter option.

Example:

urinecpt(A0 = Cl * C)  // Elimination compartment, rate proportional to concentration

NONMEM Equivalent: There isn’t a direct equivalent in NONMEM. You’d often model urine excretion using a regular compartment ($DES or built-in ADVAN subroutine) and then, if needed for steady-state calculations, you would manually ensure that the urine compartment’s contribution is handled appropriately (often by not including it in the objective function calculation at steady state).

Related Topics: deriv Statement, Compartment Models, Elimination, Steady State

Part 5: Observation and Error Models

Title: observe Statement

Keywords: observe, observation, error model, residual error, additive, proportional, combined, additiveMultiplicative, bql, censored data, M3 method, sigma

Summary: The observe statement links model predictions to observed data, defining the residual error structure. Crucially, each observe statement must include exactly one error variable. It also handles below-quantification-limit (BQL) data and provides the special sigma() function for combined error models.

Details:

Purpose: To connect model-predicted values to observed data and define how the model accounts for the difference (residual error) between them.
Syntax: observe(observedVariable = expression [, bql[ = value]][, actionCode])
- observedVariable: The name of the observed variable (e.g., CObs, EObs). This variable will be mapped to a column in your input data.
- expression: Defines the predicted value. This expression must include exactly one error variable (defined using an error statement). The form of this expression determines the error model.
- bql: (Optional) Handles below-quantification-limit (BQL) data. See “BQL Handling” for details.
- actionCode: (Optional) Allows executing sequence block code before or after the observation.
Single Error Variable Restriction: PML allows only one error variable per observe statement. NONMEM allows combining multiple EPS terms (e.g., Y = F*EXP(EPS(1)) + EPS(2)). In PML, you achieve combined error models using specific functional forms within the expression.
The sigma() Function in Combined Error Models:
- The sigma() function is a special, context-aware macro used in combined error models like additiveMultiplicative. It is not a global variable.
- Mechanism: The TDL5 translator generates separate C++ code blocks for each observe statement.
  1. Within the code block for a specific observation (e.g., CObs1), a local variable (_std) is assigned the standard deviation value from its corresponding error statement (e.g., error(CEps1=...)).
  2. The sigma() macro is simply translated to this local _std variable.
- This means that when sigma() is used inside the observe(CObs1=...) expression, it correctly and automatically refers to the standard deviation of CEps1. When used inside the observe(CObs2=...) expression, it refers to the standard deviation of CEps2. This is why you cannot use names like sigma1() or sigma2().
Common Error Models (and how to express them in PML):
- Additive: observe(CObs = C + CEps)
- Proportional: observe(CObs = C * (1 + CEps)) or observe(CObs = C * exp(CEps))
- Combined Additive and Proportional (additiveMultiplicative): This is the standard built-in form and is highly recommended. observe(CObs = C + CEps * sqrt(1 + C^2 * (CMultStdev / sigma())^2)) Here, CEps represents the additive error standard deviation, CMultStdev is an estimable parameter for the proportional standard deviation, and sigma() correctly refers to the value provided in the error(CEps=...) statement.
Multiple observe Statements: Use separate observe statements for each observed variable (e.g., one for plasma concentration, one for a PD response). Each observe statement defines its own observed variable, predicted value, and error structure.

Example (Proportional Error):

error(CEps = 0.1)  // Define the error variable and its SD
observe(CObs = C * (1 + CEps)) // Proportional error

Example (Combined Error - additiveMultiplicative):

fixef(CMultStdev = c(, 0.2, )) // Define proportional SD as a fixed effect
stparm(MultStdev = tvCMultStdev)

error(CEps = 0.05) // Define additive SD
// sigma() here will resolve to 0.05 (the value from CEps)
observe(CObs = C + CEps * sqrt(1 + C^2*(MultStdev/sigma())^2))

Example (Multiple Observations - PK and PD):

// PK observation with combined error
fixef(PK_MultStdev = c(, 0.2, ))
stparm(sPK_Mult = tvPK_MultStdev)
error(PK_Eps = 0.05)
observe(CObs = C + PK_Eps * sqrt(1 + C^2*(sPK_Mult/sigma())^2))

// PD observation with simple additive error
error(PD_Eps = 2)
observe(EObs = E + PD_Eps)

NONMEM Translation Note: When translating from NONMEM, remember that PML cannot directly combine multiple EPS terms in a single Y = ... line. You must use PML’s additiveMultiplicative form. NONMEM’s Y = F*EXP(EPS(1)) + EPS(2) is approximated in PML by this form, where sigma() provides the necessary context-specific additive standard deviation.

Related Topics: error Statement, Error Models, BQL Handling, LL Statement, Data Mapping, additiveMultiplicative

Title: `error` Statement

Keywords: error, residual error, error variable, epsilon, standard deviation, freeze, sigma

Summary: The error statement defines a residual error variable and its standard deviation. This variable is required for use within an observe statement.

Details:

Purpose: To declare a residual error variable (often with an “Eps” suffix) and specify its initial standard deviation. The standard deviation provided here is the value that the sigma() function will reference within the corresponding observe block.
Syntax: error(errorVariable[(freeze)] = standardDeviation)
- errorVariable: The name of the error variable (e.g., CEps, EEps).
- freeze: (Optional) If present, the standard deviation is fixed at the given value and not estimated.
- standardDeviation: A numeric literal (e.g., 0.1, 2.5) representing the initial estimate. Crucially, this must be a simple number. It cannot be an expression, a function call (like sqrt()), or a variable.
Default Value: If the standardDeviation is not provided, the default is 1. It is good practice to provide an explicit initial estimate.
Multiple error Statements: You need one error statement for each error variable used in your model (usually one per observe statement).
Recommended placement: Immediately before the corresponding observe statement for clarity.

Example:

error(CEps = 0.1)       // Error variable CEps, SD = 0.1
error(EEps(freeze) = 5) // Error variable EEps, fixed SD = 5

// INCORRECT - standardDeviation cannot be an expression:
// error(CEps = sqrt(0.04))

// CORRECT - Use the numeric value directly
error(CEps = 0.2) // Equivalent to sqrt(0.04)

Important Notes (CRITICAL):

Standard Deviation, NOT Variance: The value in the error statement must be the standard deviation, NOT the variance. When translating from NONMEM, you must take the square root of the $SIGMA value (if it represents variance).
- Example: NONMEM $SIGMA 0.09 corresponds to PML error(CEps = 0.3).
Numeric Literal Requirement: The standardDeviation must be a numeric literal.

NONMEM Equivalent: The error statement is conceptually similar to defining EPS terms and their corresponding variance in $SIGMA. The freeze keyword corresponds to fixing a SIGMA parameter in NONMEM. Remember the critical difference: NONMEM’s $SIGMA typically defines variance, while PML’s error statement always expects the standard deviation.

Related Topics: observe Statement, Error Models, Residual Error

Title: Error Models

Keywords: error model, residual error, additive, proportional, combined, log-additive, power, correlation, uncorrelated, LL statement

Summary: Error models describe how the model accounts for the difference between predicted and observed values. Common models include additive, proportional, combined, and power.

Details:

Purpose: To quantify the discrepancy between model predictions and observations.
Common Error Models (Review):
- Additive: observe(Obs = Pred + Eps)
- Proportional: observe(Obs = Pred * (1 + Eps))
- Log-additive: observe(CObs = C * exp(Eps))
- Combined (additive+proportional):
  - additiveMultiplicative form (Preferred): observe(CObs = C + CEps * sqrt(1 + C^2 * (CMultStdev / sigma())^2))
  - MixRatio form: observe(CObs = C + CEps * (1 + C * CMixRatio))
- Power: observe(Obs = Pred + (Pred^power)*Eps)

Correlation of Error Components and Modeling Uncorrelated Errors

The Single Epsilon Source (Perfect Correlation): In all the standard combined error models shown above, both the additive and proportional parts of the error are derived from the same single error variable (e.g., CEps). For example, in the MixRatio form, the total error term is CEps * (1 + C * CMixRatio), which expands to CEps (additive part) + CEps * C * CMixRatio (proportional part). Because both parts are driven by the exact same random draw for CEps, they are perfectly correlated.
Practical Implications: For most PK/PD applications, this perfect correlation is not a practical issue. The model’s primary goal is to describe how the total residual variance changes with the prediction, and these models accomplish that effectively.

Related Topics: observe Statement, error Statement, Residual Error, LL Statement

Title: BQL Handling (Censored Data)

Keywords: BQL, below quantification limit, censored data, M3 method, observe, bql, LLOQ, Laplacian

Summary: Below Quantification Limit (BQL) data occurs when the true value is unknown but the observed value is below the assay’s Limit of Quantification (LOQ). PML handles this using the M3 method (Laplacian estimation), calculating the likelihood that the value falls between $-\infty$ and the LOQ.

Details:

The M3 Method:
- If a sample is Quantifiable (Not BQL): The engine uses the standard Probability Density Function (PDF) based on the error model (e.g., Normal, Log-Normal).
- If a sample is Censored (BQL): The engine calculates the probability (likelihood) that the prediction plus error falls below the specified LOQ.
PML Implementation: BQL handling is activated by adding the bql option to the observe statement. There are two primary modes:
1. Dynamic BQL (Dataset-Driven):
  - Syntax: observe(Obs = Prediction + Error, bql)
  - Behavior: The censoring status and the specific LOQ value are derived from the input dataset.
  - Data Requirements: This method requires specific data preparation.
    1. BQL Flag Column: A column where 0 (or blank) indicates a valid observation, and any non-zero number (e.g., 1) indicates a censored observation.
    2. Observation Column:
      - When Flag = 0: Contains the Observed Concentration.
      - When Flag != 0: Contains the LLOQ value for that specific sample.
  - Why? The engine needs to know what the limit is to calculate the likelihood. It reads this limit directly from the observation column when the flag is raised.
2. Static BQL (Model-Driven):
  - Syntax: observe(Obs = Prediction + Error, bql = <NumericLiteral>)
  - Behavior: A fixed LOQ value is hard-coded in the model (e.g., bql = 0.5).
  - Logic: Any data value in the observation column less than or equal to this static value is automatically treated as censored.
  - Override: If a BQL flag column is explicitly mapped (see Data Mapping), the flag takes precedence over the static logic.
Mapping:
- There are no separate censor or loq mapping statements.
- BQL mapping is a sub-option of the observation mapping: obs(ModelObsVariable <- DataObsColumn, bql <- DataFlagColumn)

Example 1: Dynamic BQL (Standard Approach)

Use this when LOQ varies (e.g., different assays, different studies) or to be explicit about censoring.

PML Code:

error(CEps = 0.1)
// 'bql' keyword activates M3 method
observe(CObs = C * (1 + CEps), bql)

Column Mapping (in .cols file or R code):

// Map the Observation and the Flag
obs(CObs <- DV, bql <- BQL_Flag)

Dataset Example:

ID	Time	DV	BQL_Flag	Notes
1	1.0	12.5	0	Valid obs, conc is 12.5
1	24.0	0.5	1	Censored. 0.5 is the LLOQ

Example 2: Static BQL

Use this for quick models where the LOQ is constant and known.

PML Code:

error(CEps = 0.1)
// Any DV <= 0.05 is treated as BQL
observe(CObs = C * (1 + CEps), bql = 0.05)

Column Mapping:

// No flag mapping required (though allowed)
obs(CObs <- DV)

NONMEM Equivalent:

PML’s built-in bql option automates the IF/ELSE logic often manually coded in NONMEM for the M3 method.

$ERROR
  ; ... error definitions ...
IF (DV.GE.LLOQ) THEN
  F_FLAG=0 ; Standard PDF
  Y = CP + SD*EPS1
ELSE
  F_FLAG=1 ; M3 Likelihood
  ; PHI calculation for probability below LLOQ
  Y = PHI((LLOQ-CP)/SD)
ENDIF

Related Topics: observe Statement, Data Mapping (Column Definitions), Error Models

Title: Finding Extrema (`peak` function and alternative methods)

Keywords: peak, trough, Cmax, Tmax, Cmin, Tmin, extremum, maximum, minimum, Lagrange interpolation, peakreset, table, output, simulation

Summary: PML provides the peak function for automatically finding the maximum (peak) or minimum (trough) values of a model variable (typically concentration) and the corresponding times. An alternative, more manual method using conditional logic is also possible. Both approaches write results to output tables.

Details:

PML offers two primary methods for identifying and capturing the maximum (Cmax, Tmax) and minimum (Cmin, Tmin) values of a model variable (usually concentration, C) over time:

The peak Function (Automatic, with Interpolation):
- Purpose: Automatically finds and reports the peak (maximum) or trough (minimum) of a specified variable, along with the time at which it occurs. Uses Lagrange interpolation for more precise estimation.
- Syntax: variableName = peak(internalVariableName = expression, [max/min] [, logicalExpression])
  - variableName: A user-defined variable name to store the time of the extremum (e.g., Tmax, Tmin). This variable will be written to the output table.
  - internalVariableName: A user-defined variable to store the value. It is used internally by the peak() function and should not be declared anywhere else.
  - expression: The expression to track (usually the concentration, C).
  - max/min: (Optional) Specifies whether to find the maximum (max, default) or minimum (min).
  - logicalExpression: (Optional) A logical expression (e.g., t < 6) that restricts the search. If the expression is true, the function searches for a peak/trough as specified. If false, it searches for the opposite (trough if max, peak if min).
- Behavior:
  - Initialization: Before the first peak/trough is found (or after a peakreset), the output variables (variableName and the internalVariableName ) will be blank (missing) in the output table.
  - Detection: The peak function continuously monitors the expression and uses a 4-point window to detect potential extrema. It uses Lagrange interpolation on these points to estimate the precise time and value of the peak/trough.
  - Updating: Once a peak/trough is found, the variableName (time) and corresponding value are written to the output table. Subsequent peaks/troughs will only update these values if they are higher/lower (for max/min, respectively) than the previously found extremum.
  - If the optional logical expression is defined, the extremum is updated only when the logical expression holds true.
  - peakreset: The peakreset(internalVariableName) function resets the internal state of the peak function, causing it to start searching for a new extremum from that point forward. This is crucial for finding multiple peaks/troughs within a simulation. It is used in the sequence block.
- Restrictions:
  - Use only in the main model block, not in sequence, dobefore, or doafter blocks.
Example (Finding Cmax and Tmax):
```
Tmax = peak(Cmax = C, max) // Find the maximum concentration (Cmax) and its time (Tmax)
```
Example (Finding Cmin and Tmin within a specific time window):
```
Tmin = peak(Cmin = C, min = (t >= 44 and t <= 52)) // Find minimum between t=44 and t=52
```
Example (Multiple peaks, using peakreset):
```
sequence {
  sleep(44)  // Wait until the start of the interval of interest
  peakreset(Cmax) // Reset Cmax search
  peakreset(Cmin) // Reset Cmin search
  sleep(8)   // Observe for 8 time units
  peakreset(Cmax)
  peakreset(Cmin)
}

Tmax = peak(Cmax = C, max)
Tmin = peak(Cmin = C, min)
```
Caution: The peak function uses cubic spline interpolation, which can be sensitive to discontinuities (e.g., at the start/end of an IV infusion). Ensure sufficient simulation output density around potential extrema for accurate results.
Manual Method (Using Conditional Logic):
- Purpose: Provides more control over the extremum-finding process but requires more manual coding.
- Method: Use assignment statements and the max and min functions within the main model block, combined with the ternary operator (? :) to track the highest/lowest values and corresponding times. All variables used with this method must be initialized appropriately.
  - Initialize Cmax1 to a low value that is guaranteed to be exceeded, e.g. -1e6
  - Initialize Cmin1to a high value, e.g., 1E6.
  - Use max(current_value, previous_max) to update the maximum.
  - Use min(current_value, previous_min) to update the minimum.
  - Use the ternary operator to update the time (Tmax1, Tmin1) only when a new extremum is found.
- Important: Unlike with the peak function, you must initialize the variables used to store the maximum and minimum values.
Example (Finding Cmax and Tmax manually):
```
real(Cmax1, Tmax1, Cmin1, Tmin1)

sequence{
    Cmax1 = -1E6
    Cmin1 = 1E6
}

Cmax1 = max(C, Cmax1)
Tmax1 = (C == Cmax1 ? t : Tmax1)  // Update Tmax1 only when C equals the current Cmax1
Cmin1 = C > 0 ? min(C, Cmin1) : 1E6 // use some big value until initialization
Tmin1 = C == Cmin1 ? t : Tmin1
```
Advantages of Manual Method: Greater control, potentially more robust in some situations with discontinuities (if carefully implemented). Disadvantages of Manual Method: More code, requires careful initialization, no built-in interpolation.

Key Differences Summarized:

Feature	`peak` Function	Manual Method
Initialization	Automatic (blanks until first extremum)	Manual (must initialize Cmax1/Cmin1 to appropriate low/high values)
Interpolation	Yes (Lagrange)	No
Resetting	`peakreset(internalVariableName)`	Manual logic (typically re-initializing in a `sequence` block)
Code Complexity	Less	More
Control	Less	More
Discontinuities	Potentially sensitive (ensure sufficient output density)	Can be more robust if handled carefully
Where to Use	Main model block (not in `sequence`, `dobefore`, or `doafter`).	Main model block (variables must be declared with `real` or `double` if modified within `sequence`).
Output	Table Only	Table Only

NONMEM Equivalent: There’s no direct equivalent to PML’s peak function in NONMEM. NONMEM users typically implement custom code in $PRED or $PK to find extrema, similar to the “Manual Method” described above. NONMEM does not have built-in Lagrange interpolation for this purpose.

Title: Table Output

Keywords: table, output, results, csv, data, time points, dosepoints, covariates, observations, seq

Summary: The table statement, defined in the column definition file (or equivalent interface), specifies which variables and events should be included in the output table(s) generated by a PML model run.

Details:

Purpose: To control the generation of output tables containing simulation results, observed data, and other model-related information.
Location: The table statement is not part of the PML code itself. It’s defined in the column definition file (or equivalent interface) used by the modeling software (e.g., Phoenix NLME).
Syntax:
```
table(
  [optionalFile]
  [optionalDosepoints]
  [optionalCovariates]
  [optionalObservations]
  [optionalTimes]
  variableList
)
```
- optionalFile: Specifies the output file name (e.g., file="results.csv"). If omitted, a default file name is used.
- optionalDosepoints: Specifies that output should be generated at times when doses are administered to specific compartments (e.g., dose(A1, Aa)).
- optionalCovariates: Specifies that output should be generated when the values of specified covariates change (e.g., covr(Weight, Sex)).
- optionalObservations: Specifies that output should be generated at times when specific observations are made (e.g., obs(CObs, EObs)).
- optionalTimes: Specifies explicit time points for output generation. Can include:
  - Individual time points: time(0, 1, 2.5, 5)
  - Sequences: time(seq(0, 10, 0.1)) (generates 0, 0.1, 0.2, …, 9.9, 10)
  - Combinations: time(0, seq(1, 5, 0.5), 10)
- variableList: A comma-separated list of variables to include in the table. This can include:
  - Observed variables (CObs, EObs, etc.)
  - Predicted variables (C, E, etc.)
  - Covariates (Weight, Sex, etc.)
  - Model parameters (fixed effects, but not structural parameters directly)
  - Secondary parameters
  - User-defined variables for tracking extrema (e.g., Tmax, Cmax from the peak function, or Tmax1, Cmax1 from the manual method)
  - Time (t or mapped time variable)
  - Other derived variables
Multiple Tables: You can define multiple table statements to generate different output tables with different contents and time points.

Example:

table(file="results.csv",
      time(0, seq(1, 24, 0.5), 48),  // Output at t=0, every 0.5 from 1 to 24, and at t=48
      dose(A1),                      // Output at dose times to compartment A1
      obs(CObs),                   // Output at observation times for CObs
      covr(Weight),                // Output when Weight changes
      C, CObs, Tmax, Cmax, Tmin, Cmin  // Include these variables
     )

table(file="covariates.csv",
      time(seq(0,24, 4)),
      covr(Weight, Sex),
      Weight, Sex, Age
)

This defines two tables:

results.csv: Contains the predicted concentration (C), observed concentration (CObs), time of maximum concentration (Tmax), maximum concentration (Cmax), time of minimum concentration (Tmin) and minimum concentration (Cmin) , generated at specified times, dose times, observation times, and when Weight changes.
covariates.csv: Contains Weight, Sex and Age generated at specified times and when Weight or Sexchanges.

Important Notes:

The table statement controls output, not the model’s internal calculations.
The order of variables in the variableList determines the order of columns in the output table.
Time points specified in optionalTimes are automatically sorted.
Variables like Tmax, Cmax (from peak function or manual method) should be added in output tables, not in secondary statements.

Part 6: Delays

Title: `delay` Function

Keywords: delay, delay function, discrete delay, distributed delay, time delay, DDE, delay differential equation, gamma, Weibull, inverse Gaussian, hist

Summary: The delay function introduces a time delay into a model. It can represent either a discrete delay (all signal mediators have the same delay) or a distributed delay (delay times follow a distribution), modeling processes with a time lag.

Details:

Purpose: To model time delays, either discrete or distributed, in a dynamic system.
Syntax: delay(S, MeanDelayTime [, shape = ShapeParam][, hist = HistExpression][, dist = NameOfDistribution])
- S: The signal (expression) to be delayed. Cannot directly depend on dose-related inputs.
- MeanDelayTime: The mean delay time.
- shape = ShapeParam: (Optional) Shape parameter for the distribution (distributed delay). Interpretation depends on dist:
  - dist = Gamma or dist = Weibull: ShapeParam is the shape parameter minus 1. Must be non-negative.
  - dist = InverseGaussian: ShapeParam is the shape parameter itself.
- hist = HistExpression: (Optional) The value of S before time 0 (the “history function”). If not provided, S is assumed to be 0 before time 0.
- dist = NameOfDistribution: (Optional) The distribution: Gamma (default), Weibull, or InverseGaussian.
Discrete vs. Distributed Delay:
- If shape is not provided, a discrete delay is used: S(t - MeanDelayTime).
- If shape is provided, a distributed delay is used (convolution of S with the distribution’s PDF).
Limitations:
- Cannot be used with closed-form solutions (cfMicro, cfMacro) or matrix exponentiation.
- S cannot directly depend on dose-related inputs. Use delayInfCpt for absorption delays.
- Should be used sparingly

Example (Discrete Delay):

deriv(A = -k * delay(A, 2, hist = 10))  // A is delayed by 2 time units, initial value 10
sequence {A=10}

Example (Distributed Delay - Gamma):

deriv(A = -k * delay(A, 5, shape = 2, hist = 0, dist = Gamma)) // Gamma-distributed delay

NONMEM Equivalent: There is no direct single-function equivalent in NONMEM. Delays, especially distributed delays, often require custom coding in NONMEM using differential equations or user-defined subroutines.

Related Topics: Time Delays, Distributed Delays, gammaDelay Function, delayInfCpt Statement, Differential Equations

Title: gammaDelay Function and Transit Chain Optimization

Keywords: gammaDelay, distributed delay, gamma distribution, approximation, ODE approximation, shape parameter, optimization, performance, speed, transit compartment, ODE, systems pharmacology

Summary: Models a gamma-distributed delay using an efficient ODE approximation. It is also a powerful optimization tool for replacing manually coded transit compartment chains, often leading to significantly faster model fitting times by reducing the number of differential equations.

Details:

Purpose: The gammaDelay function has two primary uses:
1. To efficiently model a process where a time delay follows a gamma distribution.
2. To optimize complex models by replacing a computationally expensive chain of manually coded transit compartments with a single, highly optimized function.
Syntax: gammaDelay(S, MeanDelayTime, shape = ShapeParam, [, hist = HistExpression], numODE = NumberOfODEUsed)
- S: The signal (expression) to be delayed.
- MeanDelayTime: The mean of the gamma distribution (e.g., MTT).
- shape = ShapeParam: The shape parameter of the gamma distribution. When replacing a transit chain, this is the number of compartments.
- hist = HistExpression: (Optional) Value of S before time 0. Defaults to 0.
- numODE = NumberOfODEUsed: (Required) Number of ODEs for the internal approximation. Higher values are more accurate but slower. A value of 30 is often a good balance.
Accuracy and Performance: Accuracy depends on numODE.
- ShapeParam > 1: Use at least 21 ODEs.
- ShapeParam <= 1: Use at least 101 ODEs.
Limitations:
- Only for gamma-distributed delays.
- The signal to be delayed (S) cannot depend directly on dose inputs (e.g., from a dosepoint).

Advanced Usage: Optimizing Transit Compartment Chains

A chain of N identical, linear transit compartments is a common way to model a delay, and it mathematically approximates a gamma-distributed delay with shape = N. For models with many ODEs (like the CART model), these manual chains can be computationally expensive. The built-in gammaDelay() function is a more efficient way to represent this process as it reduces the number of ODEs the solver must handle.

Conversion Process:

Identify the Chain: Locate the source compartment, the destination, the number of compartments (N), and the mean transit time (MTT).
Define the “Signal”: Determine the expression for the material entering the start of the chain. This is typically (N/MTT) * source_compartment_amount.
Implement gammaDelay(): Create a new functional variable that calls gammaDelay() with the signal, where shape = N.
Refactor the deriv Statements:
- Remove the deriv statements for all intermediate transit compartments.
- Modify the source compartment’s deriv: Remove the outflow term that fed the chain.
- Modify the destination compartment’s deriv: Add the new gammaDelay variable as an input.

Example (Optimizing the CART Model):

The original CART model has a 5-compartment transit chain (A2 -> A3 -> A4 -> A5 -> A6 -> A7). We can replace this chain to speed up the model.

“Before” Code (Manual Chain):

deriv(A2  = -KTR * A2 - KC1CD * A2 + ...)
deriv(A3  = KTR * (A2 - A3))
deriv(A4  = KTR * (A3 - A4))
deriv(A5  = KTR * (A4 - A5))
deriv(A6  = KTR * (A5 - A6))
deriv(A7  = KTR * A6 + ...)
KTR = MTT/5
CTOT = A2 + A3 + A4 + A5 + A6 + A7 + ...

“After” Code (gammaDelay Implementation):

    # Signal that enters the delay chain
    signal_into_chain = (MTT/5) * A2

    # The gammaDelay function replaces the ODEs for A3, A4, A5, A6
    delayed_input_to_A7 = gammaDelay(signal_into_chain, MTT, shape = 5, numODE = 30)

    # A2's deriv no longer includes the -KTR*A2 term for the transit chain
    deriv(A2  = -KC1CD * A2 + KC1DC * A9 * F_2)

    # A7's deriv now receives the output of the gammaDelay function
    deriv(A7  = delayed_input_to_A7 + KPCEC * A7 * (1 - (CTOT / CMAXC)) * F_3 - ...)

    # A3, A4, A5, and A6 are removed from the model and from CTOT
    CTOT = A2 + A7 + A8 + A10

This reduces the number of state variables in the ODE system from 10 to 6, significantly improving performance during fitting.

Related Topics: delay Function, transit Statement, deriv Statement, Distributed Delays, Gamma Distribution, ODE Approximation

Title: transit Statement

Keywords: transit compartment, absorption, delay, absorption delay, linear chain, ntr, interpolation

Summary: Models the flow of material through a simple, linear series of linked compartments with a fixed number of stages. It is best used for the case where the number of transit is an integer, as its method for handling a non-integer number of stages is an approximation that can be inaccurate.

Details:

Purpose: transit provides a simplified way to model the flow of material through a series of linked, identical compartments. It is most suitable for describing absorption or delay processes where the number of compartments is a pre-defined integer.
Syntax: transit(dest, mtt, ntr [, in = inflow] [, out = outflow])
- dest: The final destination compartment.
- mtt: The mean transit time for the entire chain.
- ntr: A parameter representing one less than the number of transit stages. (ntr = num_stages - 1).
- in = inflow: (Optional) Additional flow into the source compartment.
- out = outflow: (Optional) Additional outflow from the dest compartment.

Handling of Integer vs. Non-Integer `ntr`

The behavior of the transit statement depends critically on whether ntr is an integer.

When ntr is an integer: The statement creates a chain of ntr + 1 compartments and is completely accurate for all inputs. This is equivalent to an Erlang distribution delay.
When ntr is a non-integer: The statement does not solve a true gamma distribution. Instead, it uses a specific logarithmic interpolation scheme between the integer numbers of compartments above and below ntr. This is an approximation.
- Disclaimer: This interpolation can lead to under-prediction of the output by up to 5% when ntr is a small non-integer (e.g., 0.5) and the dosing involves infusions or multiple boluses in close succession. This artifact makes transit a poor choice for fitting models where the number of stages is an unknown parameter to be estimated.

`transit` vs. `delayInfCpt`

This is a critical distinction for modelers:

The delayInfCpt statement is best used in situations where accuracy is important. It is based on the true gamma distribution and allows the shape parameter to be estimated robustly as a continuous, non-integer value. In addition, the delayInfCpt statement can be also used to model Weibull and inverse Gaussian absorption delay.
The transit statement is best used for Simulation: Use transit only when you have a fixed, known, integer number of compartments. In this case, transit(..., ntr = N-1, ...) is an efficient and accurate way to simulate an N-stage transit chain. Also note that the transit statement is only applicable for the gamma delay.

Important Limitations

Cannot Be Used with Closed-Form Statements: A model containing transit cannot also use cfMicro or cfMacro.

Correct Usage Example (Simulating a 3-Stage Absorption Delay)

This is the ideal use case for transit: simulating a model where it has been pre-determined that the absorption is best described by a 3-compartment chain.

test() {
    dosepoint(A1)

    # Use transit to model a 3-stage absorption process from depot to A1.
    # Since there are 3 stages, ntr = 3 - 1 = 2.
    # ntr is a fixed integer, so the result is accurate.
    transit(A1,mtt = MTT, ntr = 2, out = -Cl * C)

    C = A1 / V

    # ... observe/error statements ...
    # Parameters for a simulation, no need to fix since they are not changing:
    fixef(MTT = 2)
    fixef(Cl = 5)
    fixef(V = 50)
    # etc...
}

Related Topics: delayInfCpt Function, gammaDelay Function, Distributed Delays

Title: gammaDelay Function

Keywords: gammaDelay, distributed delay, gamma distribution, approximation, ODE approximation, shape parameter, optimization, performance, speed, transit compartment, ODE, systems pharmacology, PD delay

Summary: Models a gamma-distributed delay of an internal model signal (e.g., a PD response) by implementing an efficient Ordinary Differential Equation (ODE) approximation. It is also a powerful optimization tool for replacing manually coded transit compartment chains.

Details:

Purpose: The gammaDelay function has two primary uses:
1. To model a delayed pharmacodynamic (PD) effect or any other internal process where a time delay follows a gamma distribution.
2. To optimize complex models by replacing a computationally expensive chain of manually coded transit compartments with a single, highly optimized function.
Syntax: gammaDelay(S, MeanDelayTime, shape = ShapeParam, [, hist = HistExpression], numODE = NumberOfODEUsed)
- S: The signal (expression) to be delayed (e.g., Emax * C / (EC50 + C)).
- MeanDelayTime: The mean of the gamma distribution (e.g., MTT).
- ShapeParam: The shape parameter of the gamma distribution. When replacing a transit chain, this is the number of compartments.
- hist = HistExpression: (Optional) The value of the signal S before time 0 (the “history”). Essential for non-zero baselines. Defaults to 0.
- numODE = NumberOfODEUsed: (Required) The number of ODEs for the internal approximation.

Core Concept: An Efficient ODE Approximation

Unlike delay() function that perform direct numerical integration, gammaDelay is a high-level abstraction that generates and solves a system of numODE ordinary differential equations. As detailed by Krzyzanski (2019), this system uses properties of the binomial series to accurately and efficiently approximate the result of the gamma delay convolution.

This approach is often significantly faster than direct integration of the convolution integral, especially when the delayed signal is part of a larger ODE system. The PML engine handles the complexity of this ODE system “under the hood,” allowing the modeler to use a single, simple function call.

Choosing `numODE` for Accuracy

The accuracy of the approximation depends on the number of ODEs (numODE) and the shape parameter. Based on the work by Krzyzanski (2019): * For shape >= 1: numODE = 20 is often sufficient for high accuracy. * For shape < 1: The approximation converges more slowly, and numODE > 100 may be necessary.

Application 1: Modeling Delayed PD Responses

A primary use case is modeling the time lag between plasma concentration and the observed effect.

“Before” Code (Manual Chain):

test() {
    double(TransitionRate, Scale, lgamShape, logScale, EpsVal)
    sequence{
        TransitionRate = Shape/MeanDelayTime
        Scale = 1/TransitionRate
        lgamShape = lgamm(Shape)
        logScale = log(Scale)
        EpsVal = 1e-8
        }

    Signal = t

    deriv(x0 = Signal - TransitionRate * x0)
    deriv(x1 = Signal - (TransitionRate + 1/(t + EpsVal)) * x1)
    deriv(x2 = Signal - (TransitionRate + 2/(t + EpsVal)) * x2)
    <...> # sequential ODEs are omitted for brevity
    deriv(x199 = Signal - (TransitionRate + 199/(t + EpsVal)) * x199)
    deriv(x200 = Signal - (TransitionRate + 200/(t + EpsVal)) * x200)

    # sum
    BinomialCoeff = 1
    sum_all_x = x0
    BinomialCoeff = -BinomialCoeff * (Shape - 1)/1
    sum_all_x = sum_all_x + BinomialCoeff * x1
    BinomialCoeff = -BinomialCoeff * (Shape - 2)/2
    sum_all_x = sum_all_x + BinomialCoeff * x2
    BinomialCoeff = -BinomialCoeff * (Shape - 3)/3
    <...> # sequential sums are omitted for brevity
    BinomialCoeff = -BinomialCoeff * (Shape - 199)/199
    sum_all_x = sum_all_x + BinomialCoeff * x199
    BinomialCoeff = -BinomialCoeff * (Shape - 200)/200
    sum_all_x = sum_all_x + BinomialCoeff * x200


    # int_{0}^{t}g(t-tau)*Signal(tau)dtau
    delayedSignal_W = TransitionRate/exp(lgamShape) * (TransitionRate * t)^(Shape - 1) * sum_all_x

    stparm(MeanDelayTime = tvMeanDelayTime) # mean delay time
    stparm(Shape = tvShape) # shape parameter

    fixef(tvMeanDelayTime = c(, 2, ))
    fixef(tvShape = c(, 0.6, ))
}

“After” Code (gammaDelay Implementation):

test() {
    ##=========================================================================================================
    ##                      Define constants
    ##=========================================================================================================
    double(TransitionRate, Scale, lgamShape, logScale)
    sequence{

        # transition rate
        TransitionRate = Shape/MeanDelayTime

        # Scale
        Scale = 1/TransitionRate

        #logarithm of gamma(Shape)
        lgamShape = lgamm(Shape)

        # logaritm of Scale
        logScale = log(Scale)

        }

    #===============================================================================================
    # Model
    #===============================================================================================
    # signal to be delayed is t for t>=0, and is 0 otherwise
    Signal = t

    # using delay to find the delayed signal
    delayedSignal_gammaDelay = gammaDelay(Signal, MeanDelayTime, shape = Shape, numODE = 201)

    # analytic solution for the delayed signal
    #deriv(MeanFun = t * gammPDF(t, Shape, Scale, logScale, lgamShape))
    # Hard coded t * gammPDF(t, Shape, Scale, logScale, lgamShape) to avoid numerical difficulty
    deriv(MeanFun = TransitionRate^Shape * t^Shape / exp(lgamm(Shape)) * exp(-TransitionRate * t))
    delayedSignal_True = t * gammCDF(t, Shape, Scale, logScale, lgamShape) - MeanFun

    #===============================================================================================
    # Model parameters
    #===============================================================================================
    # structure model parameters
    stparm(MeanDelayTime = tvMeanDelayTime) # mean delay time
    stparm(Shape = tvShape) # shape parameter

    # initial value for fixed effects
    fixef(tvMeanDelayTime = c(, 2, ))
    fixef(tvShape = c(, 0.6, ))

}
}

`gammaDelay` vs. Other Delay Functions

gammaDelay vs. delayInfCpt:
- Use gammaDelay for internal model signals (e.g., PD effects).
- Use delayInfCpt for external dosing inputs from a dataset.
gammaDelay vs. transit:
- Use gammaDelay for internal model signals (e.g., PD effects).
- Use transit for external dosing inputs from a dataset.
gammaDelay vs. delay:
- gammaDelay is a specific, highly optimized function for the gamma distribution using an ODE approximation.
- delay is a more general numerical integrator that can handle other distributions (Weibull, InverseGaussian) but may be less computationally efficient than gammaDelay for the gamma case.

Related Topics: delayInfCpt Statement, transit Statement, delay Function, Distributed Delays, Absorption Delay

Part 7: Built-In Functions

Title: Built-In Mathematical Functions

Keywords: mathematical functions, sqrt, exp, log, log10, pow, abs, min, max, mod, sin, cos, tan, asin, acos, atan, sinh, cosh, tanh, asinh, acosh, atanh

Summary: PML provides a wide range of standard mathematical functions that can be used in expressions within the model. These functions operate on double-precision floating-point numbers.

Details:

Common Functions:
- sqrt(x): Square root of x.
- exp(x): Exponential function (e^x).
- log(x): Natural logarithm (base e) of x.
- log10(x): Base-10 logarithm of x.
- pow(x, y): x raised to the power of y (x^y).
- abs(x): Absolute value of x.
- min(x, y): Minimum of x and y.
- max(x, y): Maximum of x and y.
- mod(x, y): Remainder of x divided by y (modulo operation).
- sin(x), cos(x), tan(x): Trigonometric functions (input in radians).
- asin(x), acos(x), atan(x): Inverse trigonometric functions (output in radians).
- sinh(x), cosh(x), tanh(x): Hyperbolic functions.
- asinh(x), acosh(x), atanh(x): Inverse hyperbolic functions.

Example:

stparm(
  Cl = tvCl * exp(nCl),
  V  = tvV * exp(nV),
  Ka = tvKa,
  F  = ilogit(tvF)    // Example of using ilogit
)
C = A / V
Rate = Cl * C
Halflife = log(2) / (Cl / V)  // Calculate half-life using log and division

Related Topics: Expressions, stparm Statement, deriv Statement

Title: Link and Inverse Link Functions

Keywords: link function, inverse link function, ilogit, logit, probit, iprobit, iloglog, icloglog, logistic regression, categorical data

Summary: Link and inverse link functions transform variables, often to constrain them to a specific range (e.g., probabilities between 0 and 1). They are important for modeling categorical/ordinal data.

Details:

Purpose: To relate a linear predictor (which can be any real number) to a parameter that has constraints (e.g., a probability).
Common Functions:
- ilogit(x): Inverse logit function (sigmoid function). Calculates exp(x) / (1 + exp(x)). Transforms any real number x to a value between 0 and 1 (a probability).
- logit: Not directly available as a named function, use log(p/(1-p)).
- probit(p): The probit function.
- iprobit(x): Inverse probit function. Equivalent to phi(x) (CDF of the standard normal distribution).
- iloglog(x): Inverse log-log link function.
- icloglog(x): Inverse complementary log-log link function.

Example (Logistic Regression):

stparm(p = ilogit(A * time + B))  // 'p' is a probability between 0 and 1
LL(Resp, Resp * log(p) + (1 - Resp) * log(1 - p))

Related Topics: multi Statement, LL Statement, Categorical Data, Ordinal Data, Logistic Regression

Title: `if` and Ternary Operator

Keywords: if, else, ternary operator, conditional logic, sequence

Summary: Conditional logic allows different calculations or actions based on variable values. PML provides if/else within sequence blocks and the ternary operator (? :) for expressions.

Details:

if/else (within sequence blocks ONLY):
- Syntax:

sequence {
  if (condition) {
    // Statements if condition is true
  } else {
  // Statements if condition is false
  }
}

-   `condition`: A logical expression (true or false). Use C++-style
    logical operators (`&&`, `||`, `!`) within the condition.
-   *Only* usable within `sequence` blocks.

Restrictions
- if/else could be used inside sequence blocks only
- The observed variable in the LL statement cannot be used anywhere in the model outside of this statement.
Ternary Operator (for expressions):
- Syntax: condition ? value_if_true : value_if_false
- condition: A logical expression. Use nested ternary operators for complex conditions (cannot use and, or, not).
- value_if_true: Value if condition is true.
- value_if_false: Value if condition is false.
- Usable anywhere an expression is allowed (e.g., stparm, deriv, observe, LL).
  
  This is the only way to express conditional logic outside of a sequence block.
- For complex conditions (more than one else), nest ternary operators. You cannot use and, or, and not keywords. Use C++ style operators: &&, || and !.

Example:

  // CORRECT - Ternary operator outside of sequence block
  EHC = SINE <= 0 ? 0 : SINE

  // INCORRECT - if/else outside of sequence block
  EHC = if (SINE <= 0) 0 else SINE

  sequence {
  // CORRECT - if/else inside sequence block
  if (Time > 10)
  {
    DoseRate = 0
    }
  }

Example (Ternary Operator):

stparm(
  Cl = tvCl * exp(dClSex * (Sex == 1 ? 1 : 0) + nCl)  // Effect of Sex on Cl
)

This is equivalent to: If Sex equal to 1 then Cl = tvCl * exp(nClPeriod2 + nCl) otherwise Cl = tvCl * exp(nCl)

Related Topics: sequence Block, Expressions, Logical Operators

Title: Logical Operators

Keywords: logical operators, comparison operators, ==, !=, >, <, >=, <=, and, or, not

Summary: Logical operators create logical expressions (conditions) that evaluate to true or false. The allowed operators differ depending on whether they are used inside or outside a sequence block.

Details:

Comparison Operators:
- ==: Equal to
- !=: Not equal to
- >: Greater than
- <: Less than
- >=: Greater than or equal to
- <=: Less than or equal to
Logical Operators:
- Within sequence blocks: Use C++-style logical operators: && (and), || (or), ! (not).
- **Outside sequence blocks:* Nested ternary operators (? :) for conditional logic are preferred, but C++-style logical operators are also permitted.

Example:

sequence {
  if (Time > 10 && DoseRate > 0) {  // Correct
    DoseRate = 0
  }
}

stparm(
  K21 = t < 6 ? 0 : t < 8 ? tvK21 : 0  // Correct: Nested ternary
)

stparm(
  Cl = tvCl * exp(dClSex * (Sex == 1) + nCl)  // Using == for comparison
)

stparm(
  K21 = t >= 6 and t <= 8 ? tvK21 : 0  // INCORRECT: Cannot use 'and'
)

Related Topics: if and Ternary Operator, sequence Block, Expressions

Title: `sleep` Function

Keywords: sleep, sequence, time delay, action code

Summary: The sleep function, used within a sequence block, pauses execution for a specified duration, introducing time delays in action code.

Details:

Purpose: To pause the execution of a sequence block.
Syntax: sleep(duration)
- duration: Expression for the amount of time to sleep (in model time units).
sequence Block Only: sleep can only be used within a sequence block.
Relative Time: duration is relative to when sleep is encountered, not absolute time.
Stability: sleep statement should be used to ensure the stability of the algorithms

Example:

sequence {
  sleep(5)     // Pause for 5 time units
  A1 = A1 + 10  // Add 10 to compartment A1 after the delay
}

Related Topics: sequence Block, Time Delays, Action Code

Title: Statistical Distribution Functions

Keywords: statistical distributions, probability density function, PDF, cumulative distribution function, CDF, log-likelihood, lnorm, lphi, phi, dweibull, ldweibull, pweibull, lpweibull, dinvgauss, ldinvgauss, pinvgauss, lpinvgauss

Summary: PML provides functions for calculating the PDF, CDF, log PDF, and log CDF of several common statistical distributions, often used within the LL statement.

Details:

Normal Distribution:
- lnorm(x, std): Logarithm of the PDF of a normal distribution (mean 0, standard deviation std).
- lphi(x, std): Logarithm of the CDF of a normal distribution (mean 0, standard deviation std).
- phi(x): CDF of the standard normal distribution (mean 0, standard deviation 1).
Weibull Distribution:
- dweibull(x, shape, scale): PDF.
- ldweibull(x, shape, scale): Log PDF.
- pweibull(x, shape, scale): CDF.
- lpweibull(x, shape, scale): Log CDF.
Inverse Gaussian Distribution:
- dinvgauss(t, mean, shape): PDF.
- ldinvgauss(t, mean, shape): Log PDF.
- pinvgauss(t, mean, shape): CDF.
- lpinvgauss(t, mean, shape): Log CDF.
Poisson Distribution:
- lpois(mean, n): Logarithm of the probability mass function.
- ppois(mean, n): Probability mass function.
- rpois(lambda): random number
- unifToPoisson(mean, r): convert a uniform random number between 0 and 1 to a Poisson random number
Negative Binomial Distribution:
- lnegbin_rp(r, p, y): logarithm of the probability mass function of a negative binomial, distribution parameterized by r and p
- megnin_rp(r, p): generate a random sample from a negative binomial distribution parameterized by r and p
- lnegbin(mean, beta, power, y): logarithm of the probability mass function of a negative binomial distribution parameterized by mean, beta (=log(alpha)), and power
- pnegbin(mean, beta, power, y): probability mass function of a negative binomial distribution parameterized by mean, beta (= log(alpha)), and power
- rnegbin(mean, beta, power): generate a random sample from a negative binomial distribution parameterized by mean beta (= log(alpha)), and power

Example (using lnorm in an LL statement):

LL(Obs, lnorm(Obs - Pred, ErrorSD))  // Custom log-likelihood for normal distribution

Related Topics: Probability Distributions, LL Statement, Likelihood

Part 8: Additional Statements and Features

Title: Secondary Parameters

Keywords: secondary parameter, derived parameter, calculated parameter, secondary, fixed effects

Summary: Secondary parameters are quantities calculated from primary (structural) parameters, fixed effects, and other variables. They are not directly estimated but are derived. They cannot depend directly on structural parameters.

Details:

Purpose: To calculate and report derived quantities of interest (e.g., half-life, AUC, Cmax), after the main model fitting is complete. These are not part of the core model dynamics.
Syntax: secondary(parameterName = expression)
- parameterName: The name of the secondary parameter.
- expression: Defines how to calculate parameterName. Can include:
  - Fixed effects (e.g., tvV, tvCl).
  - Other secondary parameters (defined before this one).
  - Mathematical operators and functions.
  - Covariates
  - idosevar, infdosevar, infratevar variables
- Cannot Include: Structural parameters (parameters defined with stparm), random effects, or variables defined by top-level assignment. Secondary parameters are functions of fixed effects, and other derived quantities, not the individual-specific parameter values or dynamic variables.
Calculation: Calculated once after all model fitting is completed.
Multiple secondary statements: can be included.
Important Note: There is generally no direct equivalent to secondary parameters in standard NONMEM code. The secondary statement in PML is a convenience for post-processing and reporting, not for defining relationships within the core model. Do not translate simple assignments from NONMEM’s $PK or $DES blocks into secondary statements.

Example:

stparm(Cl = tvCl * exp(nCl))
stparm(V  = tvV * exp(nV))
fixef(tvCl = c(,5,))
fixef(tvV = c(,50,))

secondary(
  Halflife = log(2) / (tvCl / tvV)  // Calculate half-life from fixed effects
)

Related Topics: Parameters, stparm Statement, Fixed Effects, Calculated Values

Title: Data Mapping (Column Definitions)

Keywords: data mapping, column mapping, input data, dataset, id, time, obs, amt, rate, covariate, mdv, bql, censoring, flag, dose, addl, ss, reset, date, covr, fcovr, table

Summary: Data mapping links columns in the input dataset to variables within the PML model. This is defined in a Column Mapping file (.txt or .cols) or the data mapping interface. It handles subject identification, dosing, observations (including BQL handling), covariates, and output table configuration.

Details:

Purpose: To instruct the execution engine how to interpret specific columns in the input dataset.
Location: These statements exist outside the main model(){...} block.
Key Mapping Statements:
- id(columnName): Maps a column to the subject identifier.
  - Example: id(SubjectID)
- time(columnName): Maps a column to the time variable t. Values can be decimal numbers or hh:mm[:sec] format (e.g., 12:30, 1:30PM).
  - Example: time(Time)
- obs(modelVariable <- dataColumn [, bql <- bqlFlagColumn]): Maps a column to an observed variable defined in an observe statement.
  - Standard Mapping: obs(CObs <- Conc)
  - BQL/Censored Mapping: To handle Below Quantification Limit (BQL) data using the M3 method (Laplacian), you map the BQL flag as a sub-option.
    - Syntax: obs(CObs <- DV, bql <- BQL_Flag)
    - Requirement: The model must contain observe(CObs = ..., bql).
    - Data Structure for BQL:
      1. bqlFlagColumn: Contains 0 (or blank) for valid observations, and any non-zero numeric value (e.g., 1) for censored/BQL observations.
      2. dataColumn (The Observation):
        
        If Flag is 0: Contains the actual observed value.
        
        If Flag is non-zero: Contains the LLOQ value (Limit of Quantification) for that specific record. The engine uses this value to calculate the likelihood (probability of being between -infinity and LLOQ).
    - Note: The BQL flag mapping is optional. If not provided, all data in the observation column is treated as valid (non-censored) observations using the probability density function.
- mdv(columnName): Maps a column to the “Missing Data Value” flag.
  - Behavior: If the value is 0 or ., the observation is present. Any other value indicates the observation is missing and should be ignored.
  - Example: mdv(MDV)
- dose(dosepointName <- amountColumn [, rateColumn]): Maps a column to a dosepoint.
  - Bolus: dose(A1 <- Amt)
  - Infusion: dose(A1 <- Amt, Rate). If a rate column is mapped, it takes precedence over any rate or duration parameters defined in the PML dosepoint statement.
  - Synonyms: dose corresponds to dosepoint (or dosepoint1), dose2 corresponds to dosepoint2.
- covr(covariateName <- columnName): Maps a covariate.
  - Behavior: If a value changes at T2, the new value holds for the interval (T1, T2] (backward propagation).
  - Example: covr(Wt <- BW)
- fcovr(covariateName <- columnName): Maps a covariate with forward propagation.
  - Behavior: If a value changes at T2, the old value holds for the interval (T1, T2).
  - Note: For interpolated covariates (interpolate), covr and fcovr behave identically.
- reset(columnName=c(low, high)): Maps a reset column. If the value falls between low and high (inclusive), time and all compartments are reset to initial values.
  - Example: reset(RESET=c(1,1))
- date(columnName): Maps a date column (default format m-d-y).
- ss(columnName, doseCycleDescription): Maps a column indicating Steady State.
  - Behavior: If the column is non-zero, the system runs the specified dose cycle until steady state is reached before the current row’s event.
  - Description Syntax: Uses Reverse Polish Notation (RPN).
  - Example: ss(SS, 10 bolus(A1) 12 dt) (Give 10 to A1, wait 12 units, repeat until SS).
- addl(columnName, doseCycleDescription): Maps a column for Additional Doses.
  - Behavior: If non-zero, additional cycles are administered after the current dose.
  - Example: addl(ADDL, 24 dt 10 bolus(A1)) (Wait 24 units, then give 10 to A1).
- table(options): Defines an output table generation.
  - Syntax: table(file="name.csv", time(...), dose(...), covr(...), obs(...), variableList)
  - Function: Generates an output file with model variable values at specified times, or when specific events (doses, observations, covariate changes) occur.

Related Topics: observe Statement, BQL Handling, Covariates, dosepoint Statement, Table Output

Title: `sequence` Block

Keywords: sequence, action code, sequential execution, sleep, if, else, while, initialization, time-dependent actions

Summary: The sequence block allows for sequential execution of statements, unlike the generally declarative nature of PML. Used for initialization, time-dependent actions, and discontinuous events.

Details:

Purpose: Define statements executed in order, at specific times.
Syntax: sequence { statements }
Key Feature: Sequential Execution: Order of statements matters within a sequence block.
Statements Allowed Within sequence:
- Assignment statements (for variables modifiable when the model is stopped – integrator variables, real/double variables).
- if (test-expression) statement-or-block [else statement-or-block] (Conditional execution). Use ternary operator (? :) within expressions, not if/else.
- while (test-expression) statement-or-block (Looping)
- sleep(duration-expression) (Pauses execution)
- Function calls
Execution Timing:
- sequence blocks start before model simulation (at time 0).
- Continue until sleep or end of block.
- sleep pauses, resuming at the future time.
Multiple sequence statements can be in a model, and they are executed as if in parallel.
Reset: sequence statement(s) are restarted when a reset is encountered in the data
Restrictions:
- It is not strongly not recommended to modify inside sequence fixed effects, random effects, residual erros, structural parameters. It is impossible to modify observable variables. Example (Initializing Compartments):

sequence {
  A1 = 100  // Initialize A1 to 100 at time 0
  A2 = 0    // Initialize A2 to 0 at time 0
}

Example (Time-Dependent Action):

sequence {
  sleep(5)        // Wait 5 time units
  DoseRate = 0   // Turn off an infusion
}

Example (Loop and Conditional):

real(i)
sequence {
  i = 0
  while (i < 10) {
    if (Concentration < Threshold) {
      DoseRate = HighRate
    } else {
      DoseRate = LowRate
    }
    sleep(1)
    i = i + 1
  }
}

Related Topics: Statements, Blocks, sleep Function, if and Ternary Operator, Time-Based Models, Action Code, real Statement

Title: `real` and `double` Statements

Keywords: real, double, variable declaration, sequence, modifiable variable

Summary: real and double are synonyms, declaring variables as double-precision floating-point numbers modifiable within sequence blocks.

Details:

Purpose: Declare variables to be modified in sequence blocks.
Syntax: real(variableName1, variableName2, ...) or double(variableName1, variableName2, ...)
Double Precision: real/double variables are double-precision.
Modifiable in sequence: Can be assigned new values within sequence blocks (unlike top-level assignments).

Example:

real(Counter, Flag)  // Declare Counter and Flag

sequence {
  Counter = 0
  while (Counter < 10) {
    Counter = Counter + 1
    sleep(1)
  }
  Flag = if (Concentration > Threshold) 1 else 0 //use ternary operator
}

Related Topics: Variables, sequence Block, Assignment Statements

Part 9: Model Generation Guidelines

Title: Model Generation Guidelines

Keywords: best practices, checklist, troubleshooting, model generation, PML, common pitfalls, validation, syntax, ordering, parameterization, structural parameters, covariates, centering, error models, BQL data preparation

Summary: This checklist provides best practices and common pitfalls for generating correct PML models. It covers statement ordering, parameter definition, error models (including BQL data preparation), covariate incorporation, and built-in functions.

Details:

Statement Ordering (Recommended Order):
- Covariate Declarations: covariate, fcovariate, interpolate.
- dosepoint (and delayInfCpt): Define dosing. Before deriv referencing the dosed compartment.
- Concentration Calculation: E.g., C = A1 / V. Before deriv using C.
- deriv Statements: Define differential equations.
- error Statement(s): Define error variables (and their numeric standard deviations).
- observe Statement(s): Link predictions to observations, specify error model.
- stparm Statements: Define structural parameters. After covariate declarations, dosepoint, and calculations they depend on.
- fixef Statements: Define fixed effects (population parameters).
- ranef Statements: Define random effects (inter-individual variability).
- secondary Statements: Define secondary (derived) parameters.
- sequence Blocks: For initialization (beginning) or time-dependent actions.
Structural Parameter Definition (stparm):
- Common Styles:
  - Product * exp(eta): parameter = tvParameter * exp(nParameter) (Log-normal - for positive-only parameters like V, Cl, Ka).
  - Sum * exp(eta): parameter = (tvParameter + covariate_effects) * exp(nParameter).
  - exp(Sum + eta): parameter = exp(tvParameter + covariate_effects + nParameter).
  - ilogit(Sum + eta): parameter = ilogit(tvParameter + covariate_effects + nParameter) (For parameters between 0 and 1, like bioavailability).
  - Sum + eta: parameter = tvParameter + covariate_effects + nParameter (Normal - for parameters that can be positive or negative, like E0, Emax).
- stparm Restrictions: Within stparm, only use:
  - Fixed effects (typically tv prefix).
  - Random effects (typically n prefix).
  - Declared covariates.
  - Mathematical operators and functions.
  - Cannot use variables defined by assignment (e.g., C = A1 / V).
Covariate Incorporation:
- Declaration: Covariates must be declared (covariate, fcovariate, or interpolate).
- Centering (Continuous Covariates): Improves stability. Syntax: (CovariateName / CenterValue) or (CovariateName - CenterValue).
- Categorical Covariates: Use (CovariateName == CategoryValue) in stparm. The first category found in the data is the reference.
Error Model Specification & BQL Data Prep:
- error Statement: Define error variables and their numeric standard deviations (not variance).
- observe Statement: Link predictions to observations.
- BQL/Censored Data Pitfall (CRITICAL):
  - When using observe(..., bql) (dynamic BQL), you cannot map a separate “LOQ” column in the mapping file.
  - Data Requirement: The LOQ value must be inside the Observation (DV) column itself for any record flagged as BQL.
  - Dataset Setup:
    - Flag Column: 0 = Observed, 1 = BQL.
    - Observation Column: If Flag=0, value is Concentration. If Flag=1, value is LLOQ.
- Standard Deviation vs. Variance: Ensure the value provided in the error statement is the standard deviation. If translating from NONMEM, take the square root of $SIGMA.
Initialization:
- deriv compartments are initialized to 0. sequence overrides this.
- Only real/double, integrator, and urinecpt variables are modifiable in sequence.
Use of Built-In Functions:
- Use functions for clarity and efficiency (e.g., exp, log, sqrt, ilogit, phi, lnorm).
Review and Validation:
- Syntax Check: Double-check syntax.
- Simulation Testing: Validate by simulating scenarios.
- Parameter Estimates: Check for reasonableness.
Identifiability:
- Be aware of identifiability. Can parameters be uniquely estimated?

Example (Guideline Summary):

test() {
  fcovariate(Weight)
  fcovariate(Sex())

  dosepoint(Aa, tlag = Tlag)
  C = A1 / V

  deriv(Aa = -Ka * Aa)
  deriv(A1 = Ka * Aa - Cl * C)

  error(CEps = 0.02)
  // Dynamic BQL: The data mapped to CObs must contain the LLOQ value
  // whenever the flag mapped to 'bql' is non-zero.
  observe(CObs = C * (1 + CEps), bql)

  stparm(
    Ka   = tvKa * exp(nKa),
    V    = tvV * (Weight / 70)^dVdWt * exp(nV),
    Cl   = tvCl * exp(dCldSex * (Sex == 1)) * exp(nCl),
    Tlag = tvTlag * exp(nTlag)
  )

  fixef(
    tvKa   = c(, 1, ),
    tvV    = c(, 50, ),
    tvCl   = c(, 5, ),
    tvTlag = c(, 0.5, ),
    dVdWt  = c(, 0.75, ),
    dCldSex = c(, 0, )
  )

  ranef(diag(nV, nCl, nKa, nTlag) = c(0.09, 0.04, 0.04, 0.01))
}

Related Topics: - Basic Model Structure and Statements - Observation and Error Models (BQL Handling) - Data Mapping - Covariates

Part 10: Modeling Complex Absorption Schemes

Keywords: absorption, parallel absorption, sequential absorption, dose splitting, multiple formulations, ZOFO, lag time, tlag, duration, bioavail, split

Summary: Demonstrates how to combine multiple dosepoint statements and their options (bioavail, tlag, duration, split) to model formulations with complex absorption kinetics. It covers two main scenarios: modeling independent, mutually exclusive formulations (e.g., IV vs. Tablet) and modeling a single formulation with multiple simultaneous or sequential release pathways.

Details:

The dosepoint statement is highly flexible. By using dosepoint, dosepoint1, and dosepoint2, you can model sophisticated drug delivery systems. The correct strategy depends on your experimental design.

Method 1: Modeling Independent Formulations (Separate Dose Streams)

This is the standard and most explicit pattern for modeling different, independent routes of administration that do not occur simultaneously, such as comparing an IV dose, a Tablet dose, and a Capsule dose in separate subjects or occasions.

Key Concept: Each formulation corresponds to a separate “dose stream” defined by a dosepoint or dosepoint2 statement. Crucially, each stream is mapped to a separate dose amount column in the input dataset. This creates independent pathways for each formulation, making the model clear and the data structure intuitive. dosepoint2 is used when a second stream targets a compartment that already has a dosepoint defined.

Example: IV vs. Tablet vs. Capsule Formulations

PML Code:

test(){
    deriv(A1 = - (Cl * C) + (Aa * Ka))
    deriv(Aa = - (Aa * Ka))
    C = A1 / V

    # Define three independent dose streams
    dosepoint(A1)                   // IV Stream
    dosepoint(Aa, bioavail=F_tab)   // Tablet Stream
    dosepoint2(Aa, bioavail=F_cap)  // Capsule Stream

    # ... parameters, error models, etc. ...
    stparm(F_tab = ilogit(tvLogitFtab))
    stparm(F_cap = ilogit(tvLogitFcap))
    # ...
}

Corresponding Column Definition:

dose(A1<-"DOSEIV")
dose(Aa<-"DOSETABLET")
dose2(Aa<-"DOSECAPSULE")

Corresponding Data (each subject receives only one formulation):

ID, TIME, DOSEIV, DOSETABLET, DOSECAPSULE
1, 0, 100, ., .
2, 0, ., 100, .
3, 0, ., ., 100

Method 2: Modeling a Single Complex Formulation (Dose Splitting)

This pattern models a single, complex formulation where one administered dose is simultaneously or sequentially partitioned into multiple absorption pathways. This is the correct use case for the split argument.

A. Splitting a Single Dose into Simultaneous Pathways (`split`)

This pattern is used for a single dose amount that is administered in two different ways simultaneously (e.g., a formulation that gives an initial bolus effect while also starting a slow-release infusion).

Key Concept: This is achieved by using dosepoint and dosepoint2 and adding the crucial split argument to the first dosepoint. This allows a single dose amount column (e.g., AMT) to be mapped to both statements. The bioavail option in each statement then determines what fraction of that single AMT value goes to each pathway.

Example (50% Bolus, 50% Infusion from one dose into the same compartment):

test(){
    deriv(A1 = - (Cl * C) + (Aa * Ka))
    deriv(Aa = - (Aa * Ka))
    C = A1 / V

    # DOSE SPLITTING:
    # 'split' is ESSENTIAL to allow dosepoint2 to also use the same "AMT" column.
    dosepoint(Aa, bioavail = 0.5, split)
    dosepoint2(Aa, bioavail = 0.5, duration = 3)

    # ... observe/error and parameter definitions ...
}

Corresponding Column Definition:

dose(Aa<-"AMT", split) // Both map to the same "AMT" column
dose2(Aa<-"AMT")

B. Sequential Absorption (Zero-Order followed by First-Order - ZOFO)

This pattern models a formulation where an initial zero-order absorption phase is immediately followed by a first-order absorption phase.

Key Concept: A single dose from one AMT column is partitioned. The sequence is created by using the same parameter (D) to define the duration of the zero-order process and the tlag of the first-order process. This elegantly ensures the second process begins precisely when the first one ends.

Example Code (Sequential ZOFO):

test() {
    cfMicro(A1, Cl / V, first = (Aa = Ka))
    C = A1 / V

    # SEQUENTIAL PATHWAYS from a single dose:
    # 1. First-Order: Fraction 'Fr' is DELAYED by 'D' time units. 'split' is used.
    dosepoint(Aa, tlag = D, bioavail = Fr, split)

    # 2. Zero-Order: Fraction '(1-Fr)' is absorbed over a DURATION of 'D' time units.
    dosepoint2(A1, bioavail = 1 - Fr, duration = D)

    # ... stparm, fixef, ranef for Cl, V, Ka, D, Fr ...
}

C. Parallel Absorption with Independent Lag Times

For maximum flexibility, you can assign independent lag times to each parallel absorption pathway originating from a single dose.

Key Concept: Use the split argument along with a separate tlag in each dosepoint statement. Each tlag can be linked to a different structural parameter.

Example (Parallel ZO/FO with Independent Lags from one dose):



## **Title:** Modeling Multiple Elimination Pathways

**Keywords:** elimination, clearance, Michaelis-Menten, mixed-order,
parallel pathways, deriv, Vmax, Km, Hill

**Summary:** Describes how to model drugs that are cleared by multiple
simultaneous mechanisms (e.g., a mix of linear and saturable pathways)
by summing their individual rate equations within the `deriv` statement.

**Details:**

Pharmacokinetic models in PML are modular. If a drug is eliminated by
more than one process, the total rate of elimination is simply the sum
of the rates of each individual pathway. This is implemented by
subtracting each elimination rate term within the `deriv` statement for
the central compartment.

**Total Elimination Rate = (Rate of Pathway 1) + (Rate of Pathway 2) +
...**

**Example: Mixed Linear and Saturable (Hill-Type) Elimination**

A common scenario is a drug cleared by both a first-order (linear)
process and a capacity-limited (saturable) process.

-   **Linear Pathway Rate:** `Cl * C`
-   **Saturable Pathway Rate (with Hill):**
    `(Vmax * C^h) / (Km^h + C^h)`

**Implementation in `deriv` Statement:** The `deriv` statement for the
central compartment (`A1`) would include the drug input minus both
elimination terms.

``` pml
deriv(
    A1 = (INPUT) - (Cl * C) - (Vmax * C^h) / (Km^h + C^h)
)

Full Example Model (from previous interaction): This code shows the principle in a complex model that combines parallel absorption with mixed-order elimination. Note the two separate subtraction terms in the deriv(A1 = ...) line.

test() {
    # Parallel absorption with independent lag times
    dosepoint(Aa, bioavail = Fr_FO, tlag = Tlag_FO)
    dosepoint(A1, bioavail = 1 - Fr_FO, duration = Dur_ZO, tlag = Tlag_ZO)

    C = A1 / V
    deriv(Aa = -Ka * Aa)
    deriv(
        # The two elimination pathways are summed here by subtracting each term
        A1 = (Ka * Aa) - ( (Cl * C) + (Vmax * C^h) / (Km^h + C^h) )
    )

    # ... observe/error statements ...

    # Structural parameters including Cl, Vmax, Km, and h
    stparm(Cl   = tvCl   * exp(nCl))       // Linear clearance
    stparm(Vmax = tvVmax * exp(nVmax))     // Saturable Vmax
    stparm(Km   = tvKm   * exp(nKm))       // Saturable Km
    stparm(h    = tvh    * exp(nh))        // Hill coefficient
    # ... other stparm, fixef, ranef statements ...
}

Related Topics: deriv Statement, Compartment Models, Michaelis-Menten Kinetics

Part 11: Command-Line Execution

Title: NLME Engine Workflow Overview

Keywords: command line, execution, workflow, translate, compile, link, TDL5, C++, executable

Summary: Executing an NLME model from the command line is a multi-step process that transforms the human-readable PML model file into a runnable executable. The typical workflow involves translation of the PML code into C++, compilation and linking of the C++ code, and finally, execution of the resulting file with the specified data and engine arguments.

Details:

The process of running a PML model from a command-line interface can be broken down into four main stages:

Translation:
- The process begins with your PML model file (e.g., model.pml).
- A translator utility, such as TDL5.exe, parses the PML file.
- It converts the PML syntax into a C++ source code file (typically named Model.cpp). This step effectively translates the high-level model definition into a format that can be compiled.
Compilation:
- A C++ compiler (e.g., g++) compiles the generated Model.cpp file.
- This creates an object file (Model.o), which is the machine code representation of your specific model structure but is not yet a complete program.
Linking:
- The linker takes the compiled model object file (Model.o) and combines it with the core NLME engine libraries (e.g., libNLME7.a, libMPI_STUB.a, etc.).
- These libraries contain all the necessary functions for numerical integration, statistical calculations, optimization, and data handling.
- The output of this stage is a standalone executable file (e.g., NLME7.exe or mpiNLME7.exe for parallel versions).
Execution:
- The final executable is run from the command line.
- At this stage, you provide arguments to the executable that specify the input data files, the desired analysis method, number of iterations, output file names, and other engine settings.

This entire workflow is often managed by a scripting file (e.g., a .ps1 or .bat script) that automates these steps for the user.

Related Topics: Command-Line Syntax and Arguments, Data Input: Primary Dataset (-d1)

Title: Command-Line Syntax and Arguments

Keywords: command line, arguments, flags, syntax, options, at-file, nlmeargs

Summary: The NLME engine is controlled via command-line arguments that specify input/output files and configure the analysis. Arguments are typically preceded by a hyphen (-) or slash (/) and can be grouped into a text file for convenience.

Details:

The basic structure for executing an NLME run from the command line is:

NLME7.exe [arguments]

Where [arguments] is a space-separated list of flags and file paths.

Flag Syntax: Flags control the engine’s behavior (e.g., estimation method, number of iterations). They can be prefixed with either a hyphen (-) or a forward slash (/). For example, -n 100 and /n 100 are equivalent.
Argument File (@file): For complex runs, it is highly recommended to place arguments into a text file (e.g., nlmeargs.txt). This file can then be referenced on the command line using the @ symbol. The engine will read the file and treat each line as if it were typed on the command line. This improves readability and makes the command reusable.
Argument Order: The relative order of some arguments, particularly the data input flags (-d1, -d2, etc.), is significant. It is best practice to follow the order shown in examples.

The arguments can be categorized as follows:

Engine Configuration: Flags that control the fitting algorithm, such as estimation method (-m), number of iterations (-n), ODE solver settings (-o), etc.
Input Data Files: Flags that provide the data for the model. The primary dataset is provided with -d1. Supplemental data (e.g., for additional dosing, overriding parameters, or sequential fitting) is provided with -d2, -d3, and -d4. Each of these will be detailed in a separate topic.
Output Data Files: Flags that specify the names of output files, such as the main results file (-out_file).

Example:

Consider the following command: NLME7.exe @nlmeargs.txt -d1 cols1.txt data1.txt -out_file out000001.txt

NLME7.exe: The executable program.
@nlmeargs.txt: An argument file containing engine settings (e.g., -m 6 -n 10).
-d1 cols1.txt data1.txt: Specifies the primary input data, consisting of the column definition file cols1.txt and the data file data1.txt.
-out_file out000001.txt: Specifies that the main fitting results should be written to out000001.txt.

Related Topics: NLME Engine Workflow Overview, Data Input: Primary Dataset (-d1), Data Input: Merging Dose Data (-d2), Data Input: Overriding Fixed Effects (-d3), Data Input: Sequential Fitting Data (-d4)

Title: Data Input: Primary Dataset (-d1)

Keywords: data input, primary data, -d1, /d1, column definition, dataset

Summary: The -d1 flag is the fundamental argument used to provide the primary dataset for an NLME model fit. It requires two file paths: a column definition file that maps data columns to model variables, and the data file itself.

Details:

Every NLME command-line execution requires a primary source of data, which is specified using the -d1 flag. This dataset typically contains the subject IDs, time records, observations, covariates, and can also include dosing information.

Syntax:

-d1 <column_definition_file> <data_file>

<column_definition_file>: A text file that tells the engine how to interpret the columns in the <data_file>. It contains mapping statements like id(), time(), obs(), covr(), etc. This file acts as the bridge between your raw data and the variables defined in your PML model.
<data_file>: A text file (commonly a comma-separated value, .csv, file) containing the subject data, with columns corresponding to the mappings in the column definition file.

Example:

Consider the command: NLME7.exe -d1 cols1.txt data1.txt ...

cols1.txt (Column Definition File): This file maps the columns in data1.txt to their respective roles. A key feature is the mapping for categorical covariates like sex. The statement covr(sex<-"sex"("male" = 0, "female" = 1)) does two things: it maps the data column named “sex” to the model’s sex covariate, and it defines that the text “male” corresponds to the internal value 0, and “female” to 1.

id("id")
time("time")
dose(A1<-"dose")
covr(sex<-"sex"("male" = 0, "female" = 1))
covr(wt<-"wt")
obs(CObs<-"conc")

data1.txt (Data File): Because of the explicit mapping for sex, the data file can contain either the numeric codes or the more readable text labels. Both formats below are valid and will be processed identically by the engine.

Option A: Using Text Labels (often preferred for clarity)

##id, time, conc, dose, sex, wt
1, 0, ., 10, female, 65.6
1, 10, 1.537753, 0, female, 65.6
...

Option B: Using Numeric Codes

##id, time, conc, dose, sex, wt
1, 0, ., 10, 1, 65.6
1, 10, 1.537753, 0, 1, 65.6
...

When the engine processes the data, if the sex column contains the text female, it uses the mapping from cols1.txt to convert this to the internal value 1 for use in the model. If the column already contains the number 1, it is used directly.

Related Topics: Command-Line Syntax and Arguments, Data Input: Merging Dose Data (-d2), Data Mapping (Column Definitions)

Title: Data Input: Merging Dose Data (-d2)

Keywords: data input, merge, dosing, -d2, /d2, events, time-based, ADDL, SS

Summary: The -d2 flag allows you to provide a secondary dataset containing only dosing information. The engine merges these dosing events with the primary data from -d1, which is useful for organizing complex dosing schedules separately from observations and covariates.

Details:

For models with complex or separate dosing regimens, the -d2 flag provides a way to manage this information in its own set of files. The NLME engine reads these files and integrates the dosing events into each subject’s timeline from the primary data.

Syntax:

-d2 <column_definition_file> <data_file>

<column_definition_file>: A column definition file (cols2.txt) that maps the columns in the secondary data file. This file should define time() and one or more dosing-related statements like dose(), sscol(), addlcol(), iicol(), etc. An id() mapping is also required if subject IDs are present.
<data_file>: A data file (data2.txt) containing the dosing schedule, with columns for subject ID, time, and dose amounts/regimens.

Merge Logic and Best Practices:

The engine merges the records from -d1 and -d2 by matching subject IDs and sorting all events by time. If a dose for a subject at a specific time is defined in both data sources, the engine will administer both doses.

To avoid unintended duplicate dosing, it is a critical best practice to remove all dose-related mappings (e.g., dose(), amt(), ss(), addl()) from your primary column definition file (cols1.txt) when you are providing that information via -d2.

Example:

This example separates patient covariates and observations (-d1) from their dosing schedule (-d2).

Command: NLME7.exe -d1 cols1.txt data1.txt -d2 cols2.txt data2.txt ...

cols1.txt (Primary Column File - NO Dosing Info):

id("id")
time("time")
covr(wt<-"wt")
obs(CObs<-"conc")

data1.txt (Primary Data File - Observations/Covariates):

##id,time,conc,wt
1,0,,65.6
1,10,1.53,65.6
1,20,1.09,65.6
...

cols2.txt (Secondary Column File - Dosing ONLY):

id("id")
time("Time")
dose(A1<-"A1")

data2.txt (Secondary Data File - Dosing Schedule):

##id,A1,Time
1,100,0
1,50,24
2,100,0
...

When the engine runs, it will process subject 1 by taking their weight and concentration data from data1.txt and their dose of 100 at time 0 and dose of 50 at time 24 from data2.txt, creating a complete and correctly ordered event history for the simulation.

Related Topics: Data Input: Primary Dataset (-d1), Data Input: Overriding Fixed Effects (-d3), Data Mapping (Column Definitions)

Title: Data Input: Overriding Fixed Effects (-d3)

Keywords: data input, override, fixed effects, thetas, -d3, /d3, initial estimates, bounds, optimization, fixef

Summary: The -d3 flag provides a powerful mechanism to override the initial estimates, lower bounds, and upper bounds of fixed-effect parameters (fixef) at runtime. This allows for testing different optimization settings without recompiling the model.

Details:

The -d3 flag is designed to give the user external control over the fixef parameters (thetas) of a compiled model. This is extremely useful for exploring the parameter space, testing different starting values to aid convergence, or applying different constraints to the estimation.

Syntax:

-d3 <column_definition_file> <data_file>

<column_definition_file>: A file (cols3.txt) that defines the structure for the override data. It typically contains param() to specify the column with the parameter name, init(), low(), and high() to specify columns for the initial estimate and bounds, and map() statements to link the names in the data file to the fixef variable names in the PML model.
<data_file>: A file (data3.txt) containing the specific override values. Each row corresponds to one fixed-effect parameter.

Key Advantage:

The primary benefit of using -d3 is efficiency. You can run the exact same compiled model executable multiple times with different starting conditions or constraints simply by editing the data3.txt file. This avoids the time-consuming process of changing the PML file and re-translating/re-compiling the model for each test.

Example:

Suppose your PML model contains the following fixef statements:

fixef(tvV = c(0, 10, 100))
fixef(tvCl = c(0.001, 0.1, ))
fixef(tvV2 = c(, 10, 1000))

You can override these values at runtime using the -d3 flag with the following files:

cols3.txt (Column Definition File):

param("Parameter")
init("Initial")
low("Lower")
high("Upper")
map("tvV" <- "tvV")
map("tvCl" <- "tvCl")
map("tvV2" <- "tvV2")

data3.txt (Data File): This file sets a new initial value and new bounds for tvV, a new initial value for tvCl (while keeping its lower bound), and a new lower bound for tvV2 (while keeping its initial value and upper bound). Any field left blank is not overridden.

##Parameter,Initial,Lower,Upper
tvV,20,5,150
tvCl,0.5,0.001,
tvV2,10,1,1000

When executed with these files, the engine will start its optimization using an initial estimate of 20 for tvV, 0.5 for tvCl, etc., instead of the values hard-coded in the PML file.

Related Topics: Data Input: Primary Dataset (-d1), fixef Statement, Command-Line Syntax and Arguments

Title: Data Input: Sequential Fitting Data (-d4)

Keywords: data input, sequential fitting, -d4, /d4, post-hoc, covariates, random effects, PK/PD

Summary: The -d4 flag is a specialized argument used to merge subject-specific data to be used as covariates. Its primary purpose is to enable sequential PK/PD modeling, where individual random effect estimates from a PK model fit are supplied as known inputs to a subsequent PD model fit.

Details:

The -d4 flag provides a mechanism to feed the results of one model fit into another, which is the cornerstone of the sequential PK/PD fitting strategy.

Syntax:

-d4 <column_definition_file> <data_file>

<column_definition_file>: A file (cols4.txt) that maps the columns in the secondary data file to model variables. For sequential fitting, this file will contain an id() mapping and covr() mappings for the variables being supplied (e.g., covr(nV<-"nV")).
<data_file>: A file (data4.txt) containing the subject-specific data, such as the post-hoc random effect (eta) estimates from a previous run.

Sequential PK/PD Fitting Workflow:

This technique is used to first establish the pharmacokinetics of a drug and then, using those fixed results, estimate the pharmacodynamic parameters.

Fit PK-Only Model: First, a PK-only model is fitted to concentration-time data. The individual post-hoc estimates for the random effects (e.g., nV, nCl for each subject) are saved from this run.
Prepare PK/PD Model: A second PML model is created for the PD data. In this model, two critical changes are made:
- The random effects from the PK model (nV, nCl) are re-declared as covariates (e.g., covariate(nV, nCl)). This tells the engine that these values will be provided as input, not estimated.
- The PK fixed effects are typically frozen using the freeze keyword in the fixef statement to prevent them from being re-estimated.
Execute PK/PD Model: The PK/PD model is executed using -d1 to supply the PD observation data (e.g., Effect) and -d4 to supply the pre-calculated PK random effects.

The engine merges the data by ID, so for each subject, it uses their specific nV and nCl values from data4.txt to generate an exact individual PK profile, which then drives the PD model. This isolates the estimation process to only the unknown PD parameters.

Example:

Command: NLME7.exe -d1 cols1.txt data1.txt -d4 cols4.txt data4.txt ...

PML Snippet from PK/PD Model (test.mdl): The random effects nV and nCl are now declared as covariates.

...
fixef(tvV(freeze) = c(, 0.99, ))
fixef(tvCl(freeze) = c(, 0.98, ))
covariate(nV, nCl)
ranef(diag(nEC50, nEmax, nKe0) = c(0.04, 0.008, 0.08))
...

cols4.txt (Column Definition for -d4):

id("ID")
covr(nV<-"nV")
covr(nCl<-"nCl")

data4.txt (Data for -d4 - Post-hoc etas from PK model):

##ID, nV, nCl
1, 0.2, 0.5
2, -0.19, -0.36
...

Here, -d1 would provide the Effect vs. Time data. The -d4 files supply the known nV and nCl values for each subject, fixing their individual PK and allowing for a robust estimation of the PD parameters like tvEC50, tvEmax, etc.

Related Topics: Data Input: Primary Dataset (-d1), covariate Statement, Command-Line Syntax and Arguments

Title: Data Input Flags Summary

Keywords: data input, flags, summary, -d1, -d2, -d3, -d4, merge, override

Summary: The NLME engine uses a set of four distinct flags (-d1, -d2, -d3, -d4) to handle all data input. Each flag serves a specific purpose, from providing the primary dataset to enabling advanced sequential fitting workflows.

Details:

The following table summarizes the function and intended use of each data input flag. Understanding their distinct roles is crucial for correctly structuring a command-line model run.

Flag	Primary Purpose	Typical Data Content	Key Use Case
-d1	Provide Primary Dataset	Subject IDs, time, observations, covariates, and sometimes simple dosing regimens.	Required for all model runs. Defines the core dataset that the model will be fitted to.
-d2	Merge Supplemental Dose Events	Subject IDs, time, and detailed dosing information (e.g., `dose`, `ss`, `addl`, `rate`).	Separating complex dosing schedules. Used to keep observation/covariate data clean from a complex regimen.
-d3	Override Fixed-Effect Parameters	Parameter names (`tvV`, `tvCl`, etc.), initial estimates, and/or lower/upper bounds.	Controlling optimization without recompiling. Used to test different starting values or constraints.
-d4	Merge Subject-Specific Covariate Data	Subject IDs and values for variables to be treated as known covariates.	Sequential PK/PD fitting. Used to input individual random effect estimates from a PK fit into a PD fit.

Execution Order Note:

The relative order of these flags on the command line is significant. The engine processes them in a specific sequence. A common and safe practice is to list them in ascending order:

NLME7.exe ... -d1 ... -d2 ... -d3 ... -d4 ... -out_file ...

Related Topics: Data Input: Primary Dataset (-d1), Data Input: Merging Dose Data (-d2), Data Input: Overriding Fixed Effects (-d3), Data Input: Sequential Fitting Data (-d4)

Title: Command-Line Environment Setup

Keywords: environment, setup, PATH, compiler, libraries, prerequisites, INSTALLDIR, PML_BIN_DIR, Linux, Windows, RHEL, Ubuntu, OpenMPI

Summary: For the NLME engine and its associated tools to function correctly from the command line on either Windows or Linux, the operating system’s environment must be properly configured. This involves using environment variables to define the location of the core engine libraries, compilers, and MPI tools, and ensuring these locations are added to the system’s PATH.

Details:

The command-line workflow relies on several external programs (like the C++ compiler) and libraries. The operating system needs to know where to find these files. This is accomplished by setting specific environment variables and ensuring they are part of the system PATH.

Key Locations and Variables

NLME Core Components (INSTALLDIR):
- This is the primary environment variable, pointing to the root installation directory of the NLME engine.
- This directory contains the model translator (TDL5.exe or TDL5), platform-independent scripts, and the subdirectories containing the compiled libraries.
Platform-Specific Libraries (PML_BIN_DIR - Linux Only):
- The NLME engine provides different sets of pre-compiled libraries for different Linux distributions to ensure compatibility.
- The PML_BIN_DIR variable is used to specify which set of libraries to use. The execution script will look for libraries in the $INSTALLDIR/$PML_BIN_DIR subdirectory.
- Example: For a Red Hat Enterprise Linux (RHEL) 8/9 system, you might set PML_BIN_DIR=RHEL. For an Ubuntu 22.04/24.04 system, you would set PML_BIN_DIR=UBUNTU. If this variable is not set, the script defaults to using the libraries in the root of $INSTALLDIR.
Compiler Toolchain:
- The engine requires a specific C++ and Fortran compiler suite.
- On Windows: This is typically a specific version of MinGW64. The location is often defined by the NLMEGCCDir64 environment variable.
- On Linux: This is the system’s native GCC/GFORTRAN compiler. The execution script assumes gcc and gfortran are already available in the system’s PATH.
MPI Libraries:
- For parallel execution, the appropriate MPI implementation is required.
- On Windows: The engine links against the Microsoft MPI (MS-MPI) library. Its location can be defined by the PhoenixMSMPIDir variable.
- On Linux: The engine is based on OpenMPI. The execution script uses mpiexec and mpif90 and expects them to be in the system’s PATH. The location can be specified with the PhoenixMPIDir64 variable.

How it Works in Practice

The provided execution scripts (execNLMECmd.ps1 for Windows, execNLMECmd.sh for Linux) are designed to automate this setup. They use these environment variables to construct the correct paths and temporarily modify the PATH for the execution context, ensuring that all necessary components can be found. Users running these scripts typically only need to ensure that the initial environment variables (INSTALLDIR, PML_BIN_DIR, etc.) are correctly set.

Related Topics: NLME Engine Workflow Overview, Command-Line Syntax and Arguments, Parallel Execution (MPI)

Title: Core Engine Configuration: Method, Iterations, and ODE Solver

Keywords: method, iterations, ODE, -m, -n, -o, algorithm, FOCE, QRPEM, IT2S-EM, FO, Naive-Pooled, Laplacian, LSODA, LSODE, Runge-Kutta

Summary: The NLME engine’s core behavior is controlled by flags specifying the estimation method (-m), the maximum number of iterations (-n), and the Ordinary Differential Equation (ODE) solver (-o). These are fundamental arguments for any model fit.

Details:

These three flags are typically placed in the nlmeargs.txt file (or a similarly named at-file) and are central to defining how the model fitting is performed.

Estimation Method (`-m`)

This flag selects the statistical algorithm used to estimate the model parameters. The choice of method depends on the model type, data characteristics (e.g., continuous vs. count data), and modeling objectives.

Method Name	`-m` Flag	Description
QRPEM	`1`	Quasi-Random Parametric Expectation-Maximization. A stochastic EM algorithm.
IT2S-EM	`2`	Iterated 2-Stage Expectation-Maximization.
FOCE-LB	`3`	First-Order Conditional Estimation with Lindstrom-Bates interaction.
FO	`4`	First-Order approximation. The simplest and fastest linearization method.
FOCE-ELS / Laplacian	`5`	First-Order Conditional Estimation, ELS method. This flag covers two related methods distinguished by `-xfocehess`.
Naive-Pooled	`6`	Fits all data as if it came from a single individual, ignoring inter-individual variability.

Note on Method 5: * When -m 5 is used with -xfocehess 1 (the default), the FOCE-ELS method is selected. * When -m 5 is used with -xfocehess 0, the Laplacian method is selected.

Maximum Iterations (`-n`)

This flag sets a limit on the number of iterations the chosen estimation method will perform. This prevents excessively long run times if the model fails to converge.

Syntax: -n <number> (e.g., -n 1000)
Range: Must be a non-negative integer, typically up to 10000.
Behavior: The algorithm stops when it meets convergence criteria or when it reaches the maximum number of iterations, whichever comes first.

ODE Solver (`-o`)

For time-based models defined with deriv statements, this flag selects the numerical algorithm used to solve the system of ordinary differential equations.

Solver Name	`-o` Flag	Description
LSODE (Numerical Jac.)	`1`	Livermore Solver for ODEs. Well-suited for stiff problems. Uses a numerical Jacobian.
LSODE (Analytical Jac.)	`2`	Same as above, but requires and uses an analytical Jacobian defined in the model, which can be faster.
Runge-Kutta / DVERK	`3`	A family of non-stiff solvers.
LSODA (Numerical Jac.)	`4`	An adaptive solver that automatically switches between stiff and non-stiff methods. Uses a numerical Jacobian.
LSODA (Analytical Jac.)	`5`	Same as LSODA, but requires and uses an analytical Jacobian defined in the model.
Matrix Exponent	`6`	A direct solution method suitable for linear ODE systems. Does not apply to non-linear systems.
DOPRI5	`7`	A robust non-stiff solver (Dormand-Prince 5).

Example (nlmeargs.txt): An argument file should contain all flags on a single line, separated by spaces.

-m 5 -n 1000 -o 4

This configuration directs the engine to use the FOCE-ELS method (-m 5), run for a maximum of 1000 iterations (-n 1000), and use the adaptive LSODA solver with a numerical Jacobian (-o 4).

Related Topics: Command-Line Syntax and Arguments, Standard Error Calculation, FOCE-ELS / Laplacian Tuning, QRPEM Tuning

Title: Standard Error Calculation

Keywords: standard error, SE, Hessian, Sandwich, Fisher Score, xstderr, sand, fscore, AutoSE

Summary: The NLME engine provides multiple methods for calculating standard errors (SE), a crucial measure of parameter uncertainty. The primary methods are Hessian, Fisher Score, and the robust Sandwich method. Their availability is dependent on the chosen estimation method (-m).

Details:

Standard errors can be calculated using one of several flags. If no SE flag is specified, the engine defaults to calculating the Hessian for most methods.

The Three Main SE Methods

Hessian (-xstderr):
- Calculates SEs based on the second derivatives of the likelihood function (the Hessian matrix).
- -xstderr 1: Use the Central Difference method (Default).
- -xstderr 2: Use the Forward Difference method.
Fisher Score (-fscore):
- Calculates SEs using the Fisher Information Matrix, which is based on the first derivatives (the gradient).
- Syntax: -fscore
Sandwich (-sand):
- Calculates a robust SE estimate that “sandwiches” a variance estimate of the gradients between two Hessian estimates. It is generally more robust to model misspecification than the Hessian alone.
- Syntax: -sand

Availability by Estimation Method

The availability of each SE calculation method is strictly dependent on the estimation method (-m) being used.

Estimation Method (`-m`)	Hessian (`-xstderr`)	Fisher Score (`-fscore`)	Sandwich (`-sand`)
QRPEM (`1`)	No	Yes	No
IT2S-EM (`2`)	No	No	No
FOCE-LB (`3`)	Yes	Yes	Yes
FO (`4`)	Yes	Yes	Yes
FOCE-ELS / Laplacian (`5`)	Yes	Yes	Yes
Naive-Pooled (`6`)	Yes	Yes	Yes

Automatic Standard Error Detection (`-AutoSE`)

For convenience, especially when the success of a given SE method is uncertain, the -AutoSE flag can be used. It attempts to find a valid SE estimate by trying the methods in a specific order of preference.

Logic of -AutoSE: 1. The engine first attempts to calculate the Sandwich SE. This requires the successful calculation of both the Hessian and the Fisher Score components. 2. If the Sandwich calculation fails, the engine then checks if the Hessian calculation alone was successful. If so, it reports the Hessian-based SEs. 3. If the Hessian also failed, it finally checks if the Fisher Score calculation was successful. If so, it reports those SEs. 4. If all methods fail, the engine gives up and reports that SE calculation was unsuccessful.

Special Cases

Naive-Pooled Fallback: If a requested SE calculation fails for the Naive-Pooled (-m 6) method, the engine will automatically attempt to calculate SEs using a fallback method based on variance inflation factors.
Disabling SE Calculation: To completely disable standard error calculation for a faster run, use the flag -xstderr 0.

Related Topics: Core Engine Configuration: Method, Iterations, and ODE Solver

Title: FOCE-ELS / Laplacian Tuning and Automatic Switching

Keywords: FOCE-ELS, Laplacian, FOCE-LB, xfocehess, xlandig, xblndig, tolerance, gradTolOuter, stepTolOuter, Gaussian, non-Gaussian, BQL, LL, multi

Summary: The FOCE-ELS and Laplacian methods are powerful estimation algorithms controlled by the -xfocehess flag. The NLME engine includes crucial logic to automatically switch to the Laplacian method for models with non-Gaussian features (like BQL or custom likelihoods), even if a different FOCE method was requested.

Details:

Distinguishing FOCE-ELS from Laplacian (`-xfocehess`)

The -m 5 estimation method encompasses both FOCE-ELS and Laplacian. The choice between them is controlled by the -xfocehess flag, which determines how the Hessian matrix is calculated.

FOCE-ELS: Selected with -m 5 -xfocehess 1 (This is the default if -xfocehess is not specified).
Laplacian: Selected with -m 5 -xfocehess 0.

Automatic Switching to Laplacian

The FO, FOCE-ELS and FOCE-LB (-m 3) methods rely on assumptions that are not met by certain model types. During model translation, the engine checks for these “non-Gaussian” features.

If any of the following features are detected in the PML model, a request for FOCE-ELS or FOCE-LB will be automatically overridden, and the engine will switch to the Laplacian (-m 5 -xfocehess 0) method:

Censored Data: The model uses the bql keyword in an observe statement.
Custom Likelihoods: The model contains LL, multi, event, or count statements.
All Residual Errors Frozen: Every error() statement in the model includes the (freeze) keyword.

This automatic switch ensures that an appropriate algorithm is used for the given model structure, preventing errors and improving the chance of a successful fit.

Key Tuning and Tolerance Parameters

For fine-grained control over the FOCE/Laplacian optimization process, several flags are available to adjust the accuracy and convergence criteria of the inner (Eta) and outer (Theta/Omega) optimization loops.

Flag	Controls	Description
-xlandig	Outer Loop Accuracy (Theta/Omega)	Sets the number of significant digits of accuracy for the main parameter optimization. Default is `7`.
-xblndig	Inner Loop Accuracy (Eta)	Sets the number of significant digits of accuracy for the individual random effect optimization. Also applies to the Naive-Pooled method. Default is `13`.
-gradTolOuter	Outer Loop Gradient Tolerance	Threshold for the gradient of the outer loop. Default is `2e-4`.
-stepTolOuter	Outer Loop Step Tolerance	Threshold for the change in the parameter vector of the outer loop. Default is `1e-4`.
-gradTolInner	Inner Loop Gradient Tolerance	Threshold for the gradient of the inner loop. Default is `1.71e-5`.
-stepTolInner	Inner Loop Step Tolerance	Threshold for the change in the parameter vector of the inner loop. Default is `7.07e-8`.
-refDeltaLagl	Change in Log-Likelihood Tolerance	A convergence criterion based on minimal changes in the log-likelihood value. Default is `1e-3`.

These advanced parameters allow users to balance between computational speed and the precision of the final parameter estimates.

Related Topics: Core Engine Configuration: Method, Iterations, and ODE Solver, Standard Error Calculation

Title: QRPEM Method Tuning

Keywords: QRPEM, -m 1, stochastic, sampling, xisample, xmcpem, ximpsampdof, xscramble, burn-in, xburnin, xnonomegaburn, MAP, xmapassist

Summary: The Quasi-Random Parametric Expectation-Maximization (QRPEM) method (-m 1) is a powerful stochastic algorithm for population model fitting. Its behavior can be finely tuned with a set of parameters that control the sampling process, convergence criteria, and optimization assistance.

Details:

The QRPEM method estimates parameters by integrating over the distribution of random effects. The following flags allow for detailed control over this process.

Sampling Control

These flags determine the number and type of samples drawn during the expectation step.

Flag	Controls	Description
-xisample	Number of Samples	The number of sample points used in the quasi-random sampling process. Default is `300`. Higher values increase accuracy but also run time.
-xmcpem	Sampling Method	Switches between Quasi-Random (default) and Monte-Carlo sampling. `-xmcpem 0`: Quasi-Random (default). `-xmcpem 1`: Monte-Carlo.
-ximpsampdof	Importance Distribution	Selects the statistical distribution used for importance sampling. This can significantly impact performance, especially for non-normal posterior distributions. See table below for options. Default is Normal (`0`).
-xscramble	Scrambling Method	Applies scrambling to the quasi-random number sequence to improve its properties. `-xscramble 0`: None. `-xscramble 1`: Owen (default). `-xscramble 2`: Faure-Tezuka.

Importance Sampling Distributions (-ximpsampdof):

Value	Distribution	Notes
`0`	Multivariate Normal	Default setting.
`1`	Multivariate Double Exp.	Laplace distribution.
`2`	Direct Sampling from Prior
`3-30`	Multivariate T-dist.	The value of the flag is the degrees of freedom.
`-2`	Mixture-2 distribution
`-3`	Mixture-3 distribution

Burn-In and MAP-Assisted Optimization

These flags control preliminary phases of the optimization that help guide the algorithm to a good parameter space before the main estimation begins.

Flag	Controls	Description
-xburnin	Number of Burn-in Iterations	Specifies the number of initial iterations to run, potentially with Omegas frozen, to stabilize the estimation. Default is `0`.
-xnonomegaburn	Freeze Omegas during Burn-in	`-xnonomegaburn 1` freezes the Omega matrix (random effect variances) during the burn-in phase specified by `-xburnin`. Default is `0` (off).
-xmapassist	MAP Assistance	Activates a Maximum A Posteriori (MAP) estimation for the etas periodically during the main QRPEM optimization to assist convergence. The value of the flag sets the periodicity. Default is `0` (off).

Convergence Control

Beyond the global -n flag, these parameters provide specific convergence checks for the QRPEM algorithm.

Flag	Controls	Description
-xpemrunall	Run All Iterations	`-xpemrunall 1` forces the engine to run all iterations specified by `-n`, ignoring other convergence criteria. Default is `0`.
-emTolType	Convergence Check Type	Specifies the criteria for convergence. `0`=Default (LL & Thetas), `1`=LL & All Params, `2`=LL only, `3`=Params only. Types `1-3` use a more robust “rollout” check. Default is `0`.
-emConvLen	Convergence Check Length	The number of consecutive iterations over which the convergence criterion (`emConvCritVal`) must be met. Default is `10`.
-emConvCritVal	Convergence Critical Value	The threshold value for the convergence check. Default is `5.0`.

Related Topics: Core Engine Configuration: Method, Iterations, and ODE Solver, Standard Error Calculation

Title: Execution Modes: Simulation, Prediction, and Bootstrap

Keywords: execution mode, simulation, prediction, bootstrap, -simn, -predn, -boot, pcseed, bstrat, vpc

Summary: Beyond parameter estimation, the NLME engine can be run in several distinct modes to perform other essential tasks like simulation, prediction, or bootstrap validation. These modes are activated by specific flags that change the engine’s primary objective.

Details:

These execution modes are mutually exclusive; you can only run the engine in one of these modes at a time. They typically use the final parameter estimates from a previous model fit.

Simulation Mode (`-sim...`)

In simulation mode, the engine generates new data based on the model structure and a given set of parameters. This is useful for understanding model behavior or for generating data for a visual predictive check (VPC).

Flag	Controls	Description
-simn	Activate Simulation Mode	`-simn <number of points>` activates simulation mode and specifies the number of `ivar` (e.g., time) points to generate for the simulation output.
-simmax	Max IVAR Value	`-simmax <value>` sets the maximum value for the independent variable (e.g., time) in the simulation.
-simvar	Y-Variable Name	`-simvar <name>` specifies the name of the dependent variable (Y) in the simulation output.
-simobs	Simulate at Observations	If this flag is present, the engine will generate simulated values at the actual `ivar` points from the input data, in addition to the generated points.
-simout	Simulation Output File	`-simout <file name>` specifies the file to write the simulation results to.

Prediction Mode (`-pred...`)

In prediction mode, the engine generates model predictions for a given dataset, which can be stratified by covariates. This is the basis for generating a visual predictive check (VPC).

Flag	Controls	Description
-predn	Activate Prediction Mode	`-predn <number of replicates>` activates prediction mode and specifies the number of replicates to generate for each subject in the input data.
-predx	Prediction X-Axis	`-predx <name>` specifies the name of the x-axis variable for the prediction output.
-pstrat	Stratification Covariate(s)	`-pstrat <covariate(s)>` specifies one or more covariates on which to stratify the predictions. This is essential for covariate-based VPCs.
-pcpi / -pcpe	Prediction Intervals	`-pcpi nn-nn-nn` and `-pcpe nn-nn-nn` define the internal and external percentile intervals to be calculated and reported (e.g., `5-50-95`).
-predpc	Prediction Correction	Activates prediction correction, which helps account for discrepancies between the model and the observed data. Proportional correction is the default.
-predvpc	Variance Prediction Correction	Activates variance-prediction correction, a more advanced form of prediction correction. Implies `-predpc`.
-predout	Prediction Output File	`-predout <file name>` specifies the file to write the prediction results to.
-pcseed	Predictive Check Seed	`-pcseed <random seed>` provides a seed for the random number generator used during the predictive check, ensuring reproducibility.
-vpcnpde	NPDE Calculation	If present, instructs the engine to calculate Normalized Prediction Distribution Errors (NPDE).

Bootstrap Mode (`-boot...`)

Bootstrap is a resampling technique used to assess parameter uncertainty and model stability. The engine re-estimates the parameters on many newly generated datasets (resampled with replacement from the original data).

Flag	Controls	Description
-boot	Activate Bootstrap Mode	`-boot <random seed>` activates bootstrap mode and provides a seed for the random number generator to ensure the resampling is reproducible. The output of a bootstrap run is typically a set of parameter estimates for each replicate.
-bstrat	Stratification Covariate(s)	`-bstrat[1,2,3] <covariate name>` specifies one or more covariates on which to stratify the resampling process. This ensures that the proportion of subjects in each covariate category is maintained in the resampled datasets.

Related Topics: Core Engine Configuration: Method, Iterations, and ODE Solver

Title: Parallel Execution (MPI)

Keywords: parallel, MPI, mpiexec, performance, multi-core, speed, nodes, shared drive

Summary: The NLME engine can leverage multi-core processors to significantly decrease run times by using the Message Passing Interface (MPI) for parallel computation. Activating MPI involves using a specific executable, setting MPI-related flags, and ensuring the environment is correctly configured for distributed processing.

Details:

For computationally intensive models, especially those with large datasets or complex algorithms like QRPEM, parallel execution is essential. The NLME engine achieves this by distributing the workload (typically calculations for different subjects) across multiple processes.

Activating MPI Mode

The key to enabling parallel processing is linking the model against the appropriate MPI library during the compilation/linking stage.

The execution script (execNLMECmd.ps1 or execNLMECmd.sh) handles this automatically. If the MPIFLAG is set to MPIYES, it links against the real MPI library (libmsmpi.a on Windows, or uses mpif90 on Linux).
This creates a parallel-capable executable, typically named mpiNLME7.exe.
If MPIFLAG is MPINO, it links against a stub library (libMPI_STUB.a), creating the standard serial executable NLME7.exe.

Running in Parallel

Once the mpiNLME7.exe executable is created, it must be launched using the MPI execution wrapper, mpiexec.

Flag	Controls	Description
-n (in `mpiexec`)	Number of Nodes (Processes)	Specifies the number of parallel processes to launch. This is the most critical argument for `mpiexec`. It should generally be set to the number of available CPU cores.
-wdir (in `mpiexec`)	Working Directory	For runs on a remote cluster or shared drive, this flag tells `mpiexec` where to find the necessary input files and executable.
-localonly (Linux)	Local Host Execution	A flag in the Linux script (`execNLMECmd.sh`) to indicate that all nodes should be run on the local machine.
SHARED_DRIVE (env var)	Shared Drive Location	An environment variable used by the execution scripts on Windows to specify the path to a shared network drive where files for a distributed run are placed. This is necessary if the run spans multiple machines.

Example Command (Simplified from Script)

A typical parallel execution command constructed by the scripts looks like this:

Windows (local execution): mpiexec -n 8 .\mpiNLME7.exe [nlme_arguments]

This command launches 8 parallel processes of mpiNLME7.exe on the local machine.

Linux (shared drive execution): mpiexec -n 16 -wdir /mnt/shared/run1 ./mpiNLME7.exe [nlme_arguments]

This command launches 16 parallel processes. The -wdir flag tells mpiexec that the context for the run (the executable and its input files) is located in the /mnt/shared/run1 directory. Each of the 16 processes will start with this as its working directory.

This distributed setup is essential for running on high-performance computing (HPC) clusters where the workload is spread across multiple physical or virtual machines that all have access to a shared file system.

Related Topics: NLME Engine Workflow Overview, Command-Line Environment Setup

Title: Additional Configuration and Diagnostic Flags

Keywords: sort, csv, opcodes, modeltext, tolerance, rtol, atol, diagnostic

Summary: In addition to the main configuration flags, the NLME engine offers several other arguments to control data pre-processing, provide diagnostic outputs, and set specific numerical tolerances for the ODE solvers.

Details:

The following flags provide additional control over the engine’s behavior and can be useful for specific data formats or for debugging model execution.

Data Handling Flags

These flags control how the engine handles the input data files.

-sort: If present, this flag instructs the engine to pre-sort the input data by subject ID and then by the time variable. This is often the default behavior but can be explicitly disabled if the data is already sorted or if sorting is not desired. Note: Sorting is automatically disabled for models containing a reset statement, as the order of records is critical in that case.
-csv: If present, this flag enforces a strict comma-separated value format for the input data. If omitted, the engine will attempt to guess the delimiter (e.g., comma, space, or tab).

Output Flags

-modeltext <model-file.pml>: This flag provides the name of the original PML model file to the engine. The engine then includes the full text of the model in its primary output file (dmp.txt), which is useful for archiving and ensuring that results are always paired with the exact model code that generated them.

Numerical Tolerance Flags

The engine has many advanced flags for controlling the numerical tolerances of the various algorithms, especially the ODE solvers. While the engineParams R wrapper provides robust defaults, these can be set manually for expert users.

-rtol <value>: Sets the relative tolerance for the ODE solver (e.g., 1e-6).
-atol <value>: Sets the absolute tolerance for the ODE solver (e.g., 1e-6).
-nmxstep <value>: Sets the maximum number of internal steps the ODE solver can take between output time points (e.g., 5000).

Part 1: Basic Model Structure and Statements

Title: Basic PML Model Structure

Title: Statements

Title: Blocks

Title: Comments

Title: Variables

Title: Predefined Variables

Title: Assignment Statements

Part 2: Parameter Declarations (stparm, fixef, ranef)

The group Statement

Common Parameterization Styles

Title: fixef Statement

Title: ranef Statement

Title: Fixed Effects

Title: Random Effects

Title: Log-Normal Distribution (in Parameter Definitions)

Title: Normal Distribution (in Parameter Definitions)

Title: Bounds (in fixef statements)

Part 3: Covariates

Title: Covariates

Title: covariate Statement

Title: fcovariate Statement

Title: interpolate Statement

Title: Continuous Covariates

Title: Categorical Covariates

Title: Occasion Covariates and Inter-Occasion Variability (IOV)

Title: Time-Varying Covariates

Part 4: Structural Model Definition

Title: deriv Statement

Title: Compartment Models

Title: cfMicro Statement

Title: cfMacro and cfMacro1 Statements

Title: Micro-Constants vs. Macro-Constants

Modeling Multiple Administration Routes

1. Modeling Independent Dose Streams (Mutually Exclusive Formulations)

2. Splitting a Single Dose Amount (split argument)

Title: urinecpt Statement

Part 5: Observation and Error Models

Title: error Statement

Title: Error Models

Correlation of Error Components and Modeling Uncorrelated Errors

Title: BQL Handling (Censored Data)

Title: Finding Extrema (peak function and alternative methods)

Title: Table Output

Part 6: Delays

Title: delay Function

Advanced Usage: Optimizing Transit Compartment Chains

Handling of Integer vs. Non-Integer ntr

transit vs. delayInfCpt

Important Limitations

Correct Usage Example (Simulating a 3-Stage Absorption Delay)

Core Concept: An Efficient ODE Approximation

Choosing numODE for Accuracy

Application 1: Modeling Delayed PD Responses

gammaDelay vs. Other Delay Functions

Part 7: Built-In Functions

Title: Built-In Mathematical Functions

Title: Link and Inverse Link Functions

Title: if and Ternary Operator

Title: Logical Operators

Title: sleep Function

Title: Statistical Distribution Functions

Part 8: Additional Statements and Features

Title: Secondary Parameters

Title: Data Mapping (Column Definitions)

Title: sequence Block

Title: real and double Statements

Part 9: Model Generation Guidelines

Title: Model Generation Guidelines

Part 10: Modeling Complex Absorption Schemes

Method 1: Modeling Independent Formulations (Separate Dose Streams)

Method 2: Modeling a Single Complex Formulation (Dose Splitting)

A. Splitting a Single Dose into Simultaneous Pathways (split)

B. Sequential Absorption (Zero-Order followed by First-Order - ZOFO)

C. Parallel Absorption with Independent Lag Times

Part 11: Command-Line Execution

Title: NLME Engine Workflow Overview

Title: Command-Line Syntax and Arguments

Title: Data Input: Primary Dataset (-d1)

The `group` Statement

Title: `fixef` Statement

Title: `ranef` Statement

Title: Bounds (in `fixef` statements)

Title: `covariate` Statement

Title: `fcovariate` Statement

Title: `interpolate` Statement

Title: `deriv` Statement

Title: `cfMicro` Statement

Title: `cfMacro` and `cfMacro1` Statements

2. Splitting a Single Dose Amount (`split` argument)

Title: `urinecpt` Statement

Title: `error` Statement

Title: Finding Extrema (`peak` function and alternative methods)

Title: `delay` Function

Handling of Integer vs. Non-Integer `ntr`

`transit` vs. `delayInfCpt`

Choosing `numODE` for Accuracy

`gammaDelay` vs. Other Delay Functions

Title: `if` and Ternary Operator

Title: `sleep` Function

Title: `sequence` Block

Title: `real` and `double` Statements

A. Splitting a Single Dose into Simultaneous Pathways (`split`)

Estimation Method (`-m`)

Maximum Iterations (`-n`)

ODE Solver (`-o`)

Automatic Standard Error Detection (`-AutoSE`)

Distinguishing FOCE-ELS from Laplacian (`-xfocehess`)

Simulation Mode (`-sim...`)

Prediction Mode (`-pred...`)

Bootstrap Mode (`-boot...`)