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

End of Part 1

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

Title: `stparm` Statement

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

Summary: The stparm statement defines structural parameters in the model, which are the core parameters describing the pharmacokinetic/pharmacodynamic processes. stparm statements can include fixed effects, random effects, and covariate effects.

Details:

Purpose: Defines structural parameters and their relationships to fixed effects, random effects, and covariates.
Syntax: stparm(parameter = expression)
- parameter: The name of the structural parameter (e.g., V, Cl, Ka, EC50).
- expression: An expression defining how the parameter is calculated. This expression can include:
  - Fixed effects (typically named with a tv prefix, e.g., tvV).
  - Random effects (typically named with an n prefix, e.g., nV).
  - Covariates.
  - Mathematical operators and functions.
  - Time-dependent logic using the ternary operator (? :).
Common Distributions:
- Log-Normal: parameter = tvParameter * exp(nParameter) (Parameter is always positive)
- Normal: parameter = tvParameter + nParameter
- Logit-Normal: parameter = ilogit(tvParameter + nParameter) (Parameter is between 0 and 1)
Multiple stparm Statements: A model can have multiple stparm statements.
Important Note: stparm defines structural parameters. These are typically the parameters you are interested in estimating. Variables defined by simple assignment in NONMEM’s $PK (before the Y= line) should not be defined using stparm in PML unless they are also associated with a fixef (and thus represent an estimable parameter – a THETA in NONMEM). If a NONMEM variable is assigned a value in $PK and is not a THETA, represent it in PML with a top-level assignment statement, not with stparm.
Execution: Structural parameter statements are executed before anything else, except sequence statements to initialize them.
Conditional Logic: Use the ternary operator (? :) for conditional logic within stparm expressions, not if/else.
Important Note: When using time-dependent logic within stparm, 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))) // Clearance (log-normal) Example with occasion covariate Period with 2 levels amd categorical covariate Sex with 2 levels

Related Topics: 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.

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

End of Part 2

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 (recommended)
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

pyDarwin Automation Note: When searching for IOV on a parameter, the token should provide two options: one with only the base IIV random effect, and one that adds the IOV random effects. The token must swap out both the stparm expression and the corresponding ranef block.

Correct Token Structure for IOV on Clearance (_nCl):

"_nCl": [
    [
        "* exp(nCl)",
        "ranef(diag(nCl) = c(1))"
    ],
    [
        "* exp(nCl + (Occasion==1)*nClOccasionx1 + (Occasion==2)*nClOccasionx2)",
        "ranef(diag(nCl) = c(1))\\n\\tranef(diag(nClOccasionx1) = c(1), same(nClOccasionx2))"
    ]
]

This correctly links the structural model change to the necessary change in the random effects structure.

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

End of Part 3

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

Details:

Purpose: To define where 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).
- dobefore = sequenceStmt: (Optional) sequence block executed before dose delivery.
- doafter = sequenceStmt: (Optional) sequence block executed after dose delivery (or infusion completion).
- split: (Optional, rarely used) For UI compatibility.
- 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.
Multiple dosepoints: It is allowed to have several dosepoint statements.
dosepoint1 and dosepoint2:
- dosepoint1 is a direct synonym for dosepoint.
- dosepoint2 has a special function: it allows defining a second dosing stream on the same compartment that already has a dosepoint or dosepoint1 defined. This is used for models where a single dose is split into multiple administration profiles (e.g., part bolus, part infusion).
- split Argument: When using dosepoint and dosepoint2 on the same compartment to split a single dose amount from your data, you must add the split argument to the first dosepoint statement. 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).
  - Without split: The system would expect dose() and dose2() to map to two different amount columns (e.g., AMT1, AMT2).
  - With split: dosepoint(Aa, bioavail=0.5, split) and dosepoint2(Aa, bioavail=0.5, ...) can both be sourced from a single AMT column.
Advanced Usage: Modeling with Multiple Dosepoints A powerful feature of PML is the ability to use multiple dosepoint statements to model complex absorption. If two or more dosepoint statements exist, a single dose amount (AMT) from the input data can be programmatically split between them using the bioavail option. The sum of the bioavailability fractions across all active pathways for a given dose should typically equal 1. This is the foundation for modeling the parallel and sequential absorption schemes used for advanced drug formulations.

See “Modeling Complex Absorption Schemes” for detailed examples.

Example (Bolus dose with lag time):

dosepoint(Aa, tlag = 0.5)  // Bolus dose to Aa with a 0.5 time unit lag

Example (Zero-order infusion with bioavailability):

dosepoint(A1, duration = 2, bioavail = 0.8)  // Infusion to A1 over 2 time units, 80% bioavailability

Example (Capturing the first bolus dose to use it in the secondary statement):

dosepoint(A1, idosevar = FirstDose)

Important Note: A dosepoint statement is REQUIRED for any compartment that receives doses. Without a dosepoint statement, even if your input dataset contains AMT values for that compartment, no dosing will occur in the model.

NONMEM Equivalent: * Bolus dose: In NONMEM, you’d typically use the AMT column in your dataset along with an appropriate EVID value (usually 1) to indicate a bolus dose. * Infusion: In NONMEM, you’d use AMT, RATE (or DUR), and an EVID value (usually 4 for infusions). * tlag: Similar to NONMEM’s ALAG * bioavail: Similar to NONMEM’s F * idosevar, infdosevar, infratevar: There aren’t direct NONMEM equivalents. You might use workarounds with additional parameters or data items.

Related Topics: Dosing, Bolus Dose, Infusion, Lag Time, Bioavailability, sequence Block, Compartment Models

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

End of Part 4

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

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.

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, Resp). 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, together with the error statement, determines the error model (additive, proportional, etc.).
- bql: (Optional) Handles below-quantification-limit (BQL) data using the M3 method. See the separate entry on “BQL Handling” for details.
  - bql: Uses dynamic LLOQ, requiring a mapped CObsBQL column (or equivalent).
  - bql = <value>: Uses a static LLOQ value. Important: This must be a numeric literal, not an expression.
- actionCode: (Optional) Allows executing code (sequence block) before or after the observation (using dobefore or doafter).
Single Error Variable Restriction: This is a key difference from NONMEM. 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, not by adding multiple error variables.
Common Error Models (and how to express them in PML):
- Additive: observe(CObs = C + CEps)
  - Observed value (CObs) equals the predicted value (C) plus an error term (CEps). The error is constant regardless of the prediction’s magnitude.
- Proportional: observe(CObs = C * (1 + CEps)) or observe(CObs = C * exp(CEps))
  - The error is proportional to the predicted value. C * (1 + CEps) is a common approximation. C * exp(CEps) is a log-additive form, ensuring positivity.
- Combined Additive and Proportional: PML provides two main ways to express this:
  - additiveMultiplicative (Built-in Form): pml observe(CObs = C + CEps * sqrt(1 + C^2 * (EMultStdev / sigma())^2)) error(CEps = ...) // Define CEps as usual This form is mathematically equivalent to a combined additive and proportional error model and is often preferred for its numerical stability. EMultStdev represents the proportional error standard deviation, and sigma() represents the total variance when C=0.
  - Using a Mixing Parameter (MixRatio): pml stparm(CMixRatio = tvCMixRatio) // Define a structural parameter fixef(tvCMixRatio = c(, 0.1, )) // And its fixed effect error(CEps = ...) // Define the error variable observe(CObs = C + CEps * (1 + C * CMixRatio)) This approach is more flexible and allows for easier interpretation of the mixing parameter.
- Power: observe(Obs = Pred + (Pred^power)*Eps)
  - The error term’s magnitude scales with the predicted value raised to the power of power.
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):

error(CEps = 0.05) // Additive error SD
observe(CObs = C + CEps * sqrt(1 + C^2*(0.2/sigma())^2))  // Combined, EMultStdev=0.2

Example (Combined Error - MixRatio):

stparm(CMixRatio = tvCMixRatio)
fixef(tvCMixRatio = c(, 0.1, ))
error(CEps = 0.05)
observe(CObs = C + CEps * (1 + C * CMixRatio)) // Combined error

Example (Multiple Observations - PK and PD):

error(CEps = 0.1)
observe(CObs = C * (1 + CEps)) // PK observation

error(EEps = 2)
observe(EObs = E + EEps)  // PD observation (additive error)

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 built-in combined error forms (additiveMultiplicative or the MixRatio approach) or define a custom likelihood using the LL statement. NONMEM’s Y = F*EXP(EPS(1)) + EPS(2) is approximated in PML by the additiveMultiplicative form, as NONMEM uses a linear approximation for EXP(EPS) when the variances are small.

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

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 named with an “Eps” suffix) and specify its initial standard deviation.
Syntax: error(errorVariable[(freeze)] = standardDeviation)
- errorVariable: The name of the error variable (e.g., CEps, EEps). This name must be used in the corresponding observe statement.
- 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 for the standard deviation. Crucially, this must be a simple number. It cannot be an expression, a function call (like sqrt()), or a variable.
Default Value: If the standard deviation is not provided, the default value is 1. But not providing it is a bad practice.
Multiple error Statements: You need one error statement for each error variable used in your model (usually one per observe statement).
Recommended placement: before observe statement.

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) - standard deviation, not variance

Important Notes (CRITICAL):

Standard Deviation, NOT Variance: The standardDeviation in the error statement must be the standard deviation, NOT the variance. If you are translating a model from NONMEM, remember that the $SIGMA block often defines variances. You must take the square root of the NONMEM variance to obtain the correct standard deviation for PML.
- Example: If NONMEM has $SIGMA 0.09, the corresponding PML would be error(CEps = 0.3) (because the square root of 0.09 is 0.3).
Numeric Literal Requirement: The standardDeviation must be a numeric literal. It cannot be an expression, a variable, or a function call.

NONMEM Equivalent: The error statement combined with the error model specification in observe is conceptually similar to defining error terms in NONMEM’s $ERROR block (or within $PRED/$PREDPP). The freeze keyword corresponds to fixing the associated SIGMA parameter in NONMEM. Important Difference: NONMEM’s $SIGMA represents variance, while PML’s error statement always expects the standard deviation. You must take the square root of the variance from NONMEM’s $SIGMA when translating to PML’s error.

Related Topics: observe Statement, Error Models, Residual Error

Title: Error Models

Keywords: error model, residual error, additive, proportional, combined, additiveMultiplicative, log-additive, power

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, reflecting measurement error, biological variability, and model misspecification.
Common Error Models (Review - with emphasis on PML implementation):
- Additive: Observed = Predicted + Error
  - PML: observe(Obs = Pred + Eps)
  - Error is constant, regardless of prediction.
- Proportional: Observed = Predicted * (1 + Error) or Observed = Predicted * exp(Error)
  - PML: observe(Obs = Pred * (1 + Eps)) (approximation) or observe(Obs = Pred * exp(Eps)) (log-additive)
  - Error is proportional to the prediction.
- Log-additive: observe(CObs = C*exp(Eps))
- Combined Additive and Proportional:
  - PML (Preferred - additiveMultiplicative): pml observe(CObs = C + CEps * sqrt(1 + C^2 * (EMultStdev / sigma())^2)) error(CEps = ...) // Define CEps as usual Where ‘EMultStdev’ is proportional error standard deviation, and ‘sigma()’ represents the variance when C=0.
  - PML (Using MixRatio): pml stparm(CMixRatio = tvCMixRatio) fixef(tvCMixRatio = c(, 0.1, )) error(CEps = ...) observe(CObs = C + CEps * (1 + C * CMixRatio))
- Power: Observed = Predicted + (Predicted^power) * Error
  - PML: observe(Obs = Pred + (Pred^power)*Eps)
Choosing an Error Model: Select a model that reflects how variability changes with the magnitude of the predicted value.

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

Title: BQL Handling

Keywords: BQL, below quantification limit, censored data, M3 method, observe, bql, CObsBQL, LOQ

Summary: Below Quantification Limit (BQL) data occurs when the true value is below the assay’s detection limit. PML handles BQL data using the M3 method within the observe statement.

Details:

Problem: BQL data are censored; we only know the value is below a limit (the LOQ).
M3 Method: Standard approach (Beal, 2001). Calculates the probability of the true value being below the LOQ.
PML Implementation: The observe statement’s bql option implements the M3 method.
Two Ways to Use bql:
- observe(Obs = ..., bql) (Dynamic LOQ):
  - Creates a derived column ObsBQL (e.g., CObsBQL).
  - Optionally map a column to this flag.
  - ObsBQL = 0 (or missing): Observation is above the LOQ.
  - ObsBQL = non-zero: Observation is at or below the LOQ. The Obs value becomes the LOQ.
  - If ObsBQL is not mapped, then the values in Obs are used.
- observe(Obs = ..., bql = <value>) (Static LOQ):
  - Defines a fixed LOQ value (numeric literal).
  - Obs values less than <value> are treated as censored.
  - Mapping ObsBQL is optional; mapped values override the static value.
Mapping:
- censor(CObsBQL)
- loq(LOQ)

Example (Dynamic LOQ):

error(CEps = 0.1)
observe(CObs = C * (1 + CEps), bql)  // Dynamic BQL
// ... and in the column mappings:
// censor(CObsBQL) // Optional, for explicit BQL flags
// loq(LOQ)

Example (Static LOQ):

error(CEps = 0.1)
observe(CObs = C * (1 + CEps), bql = 0.05)  // Static LOQ of 0.05

NONMEM Equivalent: In NONMEM, you’d use the M3 method.

$ERROR
  CP=A(2)/V
  PROP=CP*RUVCV
  ADD=RUVSD
  SD=SQRT(PROP*PROP+ADD*ADD)
IF (DV.GE.LLOQ) THEN
  F_FLAG=0 ; ELS
  Y=CP + SD*EPS1
ELSE
  F_FLAG=1 ; LIKELIHOOD
 Y=PHI((LLOQ-CP)/SD))
ENDIF

Related Topics: observe Statement, Error Models, Censored Data, Data Mapping

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.

End of Part 5

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

Keywords: gammaDelay, distributed delay, gamma distribution, approximation, ODE approximation, shape parameter

Summary: The gammaDelay function models a gamma-distributed delay using an ODE approximation, which can be faster than the delay function for complex models.

Details:

Purpose: Efficiently model a delay where the delay time follows a gamma distribution.
Syntax: gammaDelay(S, MeanDelayTime, shape = ShapeParam, [, hist = HistExpression], numODE = NumberOfODEUsed)
- S: The signal to be delayed.
- MeanDelayTime: The mean of the gamma distribution.
- shape = ShapeParam: The shape parameter of the gamma distribution. Not shape parameter minus one, unlike delay.
- hist = HistExpression: (Optional) Value of S before time 0. Defaults to 0.
- numODE = NumberOfODEUsed: (Required) Number of ODEs for the approximation. Higher values are more accurate but slower. Maximum value is 400.
ODE Approximation: gammaDelay approximates the convolution integral with a system of ODEs.
Accuracy and Performance: Accuracy depends on numODE.
- ShapeParam > 1: Use at least 21 ODEs.
- ShapeParam <= 1: Use at least 101 ODEs.
Advantages over delay: Can be significantly faster than delay(..., dist=Gamma).
Limitations:
- Only for gamma-distributed delays.
- The signal to be delayed cannot depend directly on dose inputs

Example:

deriv(A = -k * gammaDelay(A, 5, shape = 3, numODE = 30))  // Gamma-distributed delay
sequence {A=10}

NONMEM Equivalent: There’s no direct equivalent in NONMEM. You’d typically implement a gamma-distributed delay using a series of transit compartments, which approximates a gamma distribution.

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

Title: `delayInfCpt` Statement

Keywords: delayInfCpt, absorption delay, distributed delay, compartment, dosepoint, inflow, outflow, gamma, Weibull, inverse Gaussian

Summary: The delayInfCpt statement models a distributed delay specifically for input into a compartment. This is the correct way to model absorption delays with a distributed delay time.

Details:

Purpose: Model distributed delays in the absorption process (or any input into a compartment). Necessary because delay cannot handle dose-related inputs directly.
Syntax: delayInfCpt(A, MeanDelayTime, ParamRelatedToShape [, in = inflow][, out = outflow][, dist = NameOfDistribution])
- A: The compartment receiving the delayed input. Can receive doses via dosepoint.
- MeanDelayTime: Mean of the delay time distribution.
- ParamRelatedToShape: Related to the shape parameter. Depends on dist:
  - dist = InverseGaussian: ParamRelatedToShape is the shape parameter.
  - dist = Gamma or dist = Weibull: ParamRelatedToShape is the shape parameter minus 1. Must be non-negative.
- in = inflow: (Optional) Additional inflow into compartment A that should also be delayed.
- out = outflow: (Optional) Outflow from compartment A that should not be delayed.
- dist = NameOfDistribution: (Optional) Delay time distribution: Gamma (default), Weibull, or InverseGaussian.
Relationship to dosepoint: Used with a dosepoint statement for compartment A. dosepoint handles dosing, delayInfCpt models the delay.
History function: is assumed to be zero

Example (One-compartment model with gamma-distributed absorption delay):

delayInfCpt(A1, MeanDelayTime, ShapeParamMinusOne, out = -Cl * C)
dosepoint(A1)
C = A1 / V

Example (Two-compartment model, two absorption pathways, each with gamma delay):

delayInfCpt(Ac1, MeanDelayTime1, ShapeParamMinusOne1, out = -Cl * C - Cl2 * (C - C2))
dosepoint(Ac1, bioavail = frac)  // Fraction 'frac' goes through pathway 1
delayInfCpt(Ac2, MeanDelayTime2, ShapeParamMinusOne2)
dosepoint(Ac2, bioavail = 1 - frac) // Remaining fraction goes through pathway 2
deriv(A2 = Cl2 * (C - C2))
C = (Ac1 + Ac2) / V
C2 = A2 / V2

NONMEM Equivalent: There isn’t a single, direct equivalent in NONMEM. You would likely use a combination of: * An absorption compartment. * A series of transit compartments (to approximate a distributed delay, particularly a gamma distribution). * Potentially, user-defined subroutines (for more complex delay distributions).

Related Topics: Absorption Delay, Distributed Delays, delay Function, gammaDelay Function, dosepoint Statement, Compartment Models

Title: `transit` Statement Keywords: transit compartment, absorption

Summary: The transit statement models the flow of material through a series of linked compartments.

Details: * Purpose: transit statement can be used to model the flow of material through a series of linked compartments * Syntax: transit(dest, source, mtt, num [, in = inflow] [, out = outflow]) * dest: destination compartment * source: source compartment, could be a dosepoint * mtt: mean transit time * num: number of transit compartments * in = inflow: (Optional) Specifies any additional inflow into destination compartment * out = outflow: (Optional) Specifies any outflow from destination compartment

Restrictions: Cannot be used in a model with any closed-form statement Example:

transit(A1, A0, mtt, num, out = - Cl * A1/V)

NONMEM Equivalent: The transit statement is conceptually similar to setting up a series of compartments with first-order transfer between them in NONMEM, but it handles the underlying equations more implicitly.

Related Topics: Compartment Models

End of Part 6

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

End of Part 7

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, evid, rate, covariate, censor, loq, dvid, dose, addl, ii, table, ss, reset, date

Summary: Data mapping links columns in your input dataset to variables and contexts within your PML model. This is done outside the PML code, using a column definition file or interface.

Details:

Purpose: To tell the execution engine how to interpret the data in your input file.
Column Definition File: A text file (or equivalent interface) with mapping statements.
Key Mapping Statements:
- id(columnName): Maps a column to the subject ID. Example: id(SubjectID).
- time(columnName): Maps a column to the time variable (t). Example: time(Time).
- obs(observedVariable <- columnName): Maps a column to an observed variable. Example: obs(CObs <- Conc).
- amt(columnName): Maps a column to the dose amount (with dosepoint). Example: amt(DoseAmt).
  - Can map a single column containing both bolus and infusion amounts.
- rate(columnName): Maps a column to the infusion rate (with dosepoint). Example: rate(InfRate).
- covr(covariateName <- columnName): Maps a column to a covariate (covariate or fcovariate). Example: covr(Weight <- BW).
  - Categorical Covariates:
    - With labels: covr(Sex <- SexColumn("Male" = 0, "Female" = 1))
    - Without labels: covr(Sex <- SexColumn()) (First unique value is the reference category).
- fcovr(covariateName <- columnName): Same as covr, but for fcovariate.
- censor(columnName): Maps a column to the BQL flag from observe. Example: censor(CObsBQL).
- loq(columnName): Maps a column to provide lower limits of quantification. Example: loq(LLOQ).
- mvd(columnName): Maps a column to indicate missing data values for observations. Example: mvd(MDV).
- evid(columnName): Maps a column to an event identifier. Example: evid(EVID).
- addl(columnName, doseCycleDescription): Maps a column to indicate additional doses. Example: addl(ADDL, 24 dt 10 bolus(A1)).
- ii(columnName): Maps a column to the interdose interval. Often derived from addl.
- dose(doseVariable <- columnName, cmt = compartmentNumber): Conditional mapping of doses. Maps columnName to doseVariable only when cmt matches compartmentNumber. Essential for multiple dosepoints. Example: dose(AaDose <- AMT, cmt = 1).
- cmt(columnName): Maps a column specifying the compartment number (with conditional dose). Example: cmt(CMT).
- ss(ssColumnName, doseCycleDescription): Maps a column that indicates steady-state dosing.
- reset(resetColumnName=c(lowValue, highValue)): Maps a column to reset time and compartments
- date(dateColumnName[, formatString [, centuryBase]]): maps date column
- dvid(columnName): Maps a column to an observation type identifier (multiple observation types in same data column). Example: dvid(TYPE). Used with conditional observation mapping.
- table(...): (Not for input, defines output tables).
Conditional Mapping: Use conditional logic in mappings to handle rows differently (e.g., based on CMT, DVID). Map the same data column to different variables, depending on context.

Example (Conceptual Column Mappings - PK/PD, Multiple Doses, BQL):

id(ID)
time(TIME)
amt(AMT)
addl(ADDL)        // Additional doses
ii(II)            // Interdose interval
mvd(MDV)
evid(EVID)
obs(CObs <- DV)    // PK observations
obs(EObs <- DV)  // PD observations
covr(Weight <- WEIGHT)
covr(Sex <- SEX("Male" = 0, "Female" = 1))
censor(CObsBQL)  // BQL flag for PK observations
loq(LOQ)
dose(AaDose <- AMT, cmt = 1)   // Map AMT to AaDose when CMT=1 (bolus)
dose(AaInfDose <- AMT, cmt = 1) // Map AMT to AaInfDose when CMT=1 (infusion)
rate(AaInfRate <- RATE)     // Infusion rate
dvid(TYPE) // DVID for distinguishing PK/PD observations

Related Topics: Input Data, observe Statement, dosepoint Statement, Covariates, BQL Handling, Conditional Mapping

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 not recommended to modify inside sequence fixed effects, random effects, residual errors, 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

End of Part 8

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

Summary: This checklist provides best practices and common pitfalls for generating correct PML models. It covers statement ordering, parameter definition, error models, covariate incorporation, and built-in functions, helping to avoid mistakes and improve identifiability.

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).
- Choosing the Right Style:
  - Positive-Only (V, Cl, Ka, etc.): Product * exp(eta) is generally preferred.
  - Positive or Negative (E0, Emax, etc.): Sum + eta is appropriate.
  - Between 0 and 1 (Bioavailability): ilogit(Sum + eta) is ideal.
- 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 and reduces correlation.
- Syntax: (CovariateName / CenterValue) (multiplicative) or (CovariateName - CenterValue) (additive).
- Center Value: Mean, median, or clinically relevant value.
Categorical Covariates: Use (CovariateName == CategoryValue) in stparm. The first category is the reference. Mapping is done in the column definition file.

Examples (stparm Styles with Covariates):

# Product * exp(eta) (most common for positive parameters)
stparm(V = tvV * (Weight / 70)^dVdWt * exp(nV))  # Allometric scaling
stparm(Cl = tvCl * exp(dCldSex * (Sex == 1)) * exp(nCl)) # Sex effect

# Sum * exp(eta)
stparm(V = (tvV + dVdWt * (Weight - 70)) * exp(nV))

# exp(Sum + eta)
stparm(V = exp(tvV + dVdWt * (Weight - 70) + nV))

# ilogit(Sum + eta) (for parameters between 0 and 1)
stparm(F = ilogit(tvF + dFdWt * (Weight - 70) + nF))

# Sum + eta (for parameters that can be positive or negative)
stparm(E0 = tvE0 + dE0dWt * (Weight - 70) + nE0)

fcovariate(Weight)
fcovariate(Sex()) #categorical covariate

Error Model Specification:
- error Statement: Define error variables and their numeric standard deviations.
- observe Statement: Link predictions to observations, specify error model (additive, proportional, combined, power, log-additive). Includes one error variable.
- LL Statement: For custom likelihoods.
- Use prescribed error variable names (like CEps)
- Select an error model appropriate for your data.
- “Standard Deviation vs. Variance: Ensure the value provided in the error statement is the standard deviation, not the variance. If translating from NONMEM, remember to take the square root of the $SIGMA value (if it represents variance).”
Initialization:
- deriv compartments are initialized to 0. sequence overrides this.
- Only real/double, integrator, and urinecpt variables are modifiable in sequence.
- Do not include explicit compartment initializations within sequence blocks unless specifically required
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 and SEs: Check for reasonableness. Large values or high correlations can indicate problems.
- Documentation: Include comments.
Identifiability:
- Be aware of identifiability. Can parameters be uniquely estimated? Overparameterization can lead to non-identifiability.

Example (Guideline Summary):

test() {
  fcovariate(Weight)  # Declare covariate early
  fcovariate(Sex())    # Declare categorical covariate

  dosepoint(Aa, tlag = Tlag)  # Dosepoint (with structural lag time)
  C = A1 / V                # Define concentration *before* using it

  deriv(Aa = -Ka * Aa)      # Absorption compartment
  deriv(A1 = Ka * Aa - Cl * C - Cl2 * (C - C2) - Cl3 * (C - C3))  # Central
  deriv(A2 = Cl2 * (C - C2))                                    # Peripheral 2
  deriv(A3 = Cl3 * (C - C3))                                    # Peripheral 3
  C2 = A2/V2
  C3 = A3/V3

  error(CEps = 0.02)       # Numeric SD value!
  observe(CObs = C * (1 + CEps))  # Proportional error

  # Structural parameters (log-normal, with covariate effects)
  stparm(
    Ka   = tvKa * exp(nKa),
    V    = tvV * (Weight / 70)^dVdWt * exp(nV),         # Allometric scaling
    Cl   = tvCl * exp(dCldSex * (Sex == 1)) * exp(nCl),  # Sex effect
    Cl2  = tvCl2 * exp(nCl2),
    V2   = tvV2 * exp(nV2),
    Cl3  = tvCl3 * exp(nCl3),
    V3   = tvV3 * exp(nV3),
    Tlag = tvTlag * exp(nTlag)
  )

  # Fixed effects
  fixef(
    tvKa   = c(, 1, ),
    tvV    = c(, 50, ),
    tvCl   = c(, 5, ),
    tvCl2  = c(, 2, ),
    tvV2   = c(, 30, ),
    tvCl3  = c(, 1, ),
    tvV3   = c(, 20, ),
    tvTlag = c(, 0.5, ),
    dVdWt  = c(, 0.75, ),  # Fixed effect for Weight on V
    dCldSex = c(, 0, )     # Fixed effect for Sex on Cl
  )

  # Random effects
    ranef(diag(nV, nCl, nCl2, nV2, nCl3, nV3, nKa, nTlag) =
      c(0.09, 0.04, 0.04, 0.09, 0.04, 0.09, 0.04, 0.01)
  )
}

Related Topics:

Basic Model Structure and Statements
Parameter Declarations
Covariates
Structural Model Definition
Observation and Error Models
Delays
Built-In Functions
Data Mapping

End of Part 9

Part 10: Modeling Complex Absorption Schemes

Keywords: absorption, parallel absorption, sequential absorption, dose splitting, dual absorption, ZOFO, lag time, tlag, duration, bioavail

Summary: Demonstrates how to combine multiple dosepoint statements and their options (bioavail, tlag, duration) to model formulations with complex absorption kinetics, such as simultaneous or sequential processes.

Details:

The dosepoint statement is highly flexible. By using multiple dosepoint statements, you can model complex drug delivery systems where a single dose (AMT from the data) is split into different pathways or scheduled to occur sequentially.

A. Simultaneous (Parallel) Absorption (Zero-Order and First-Order)

This pattern models a formulation with both immediate and sustained release components that begin absorbing at the same time. A single dose is split between a first-order and a zero-order process using a fractional bioavailability parameter (Fr).

Key Concept: Two dosepoint statements are used. One (dosepoint(Aa, bioavail = Fr)) directs a fraction of the dose to a first-order process. The other (dosepoint(A1, bioavail = 1 - Fr, duration = Dur)) directs the remaining fraction to a zero-order process.

Example Code (Simultaneous ZO/FO):

test() {
    # PARALLEL PATHWAYS:
    # 1. First-order absorption for fraction 'Fr'.
    dosepoint(Aa, bioavail = Fr)
    # 2. Zero-order absorption for fraction '(1 - Fr)' over duration 'Dur'.
    dosepoint(A1, bioavail = 1 - Fr, duration = Dur)

    # Model structure
    C = A1 / V
    deriv(Aa = -Ka * Aa)
    deriv(A1 = Ka * Aa - Cl * C)
    
    # ... observe/error statements ...

    # Structural parameters
    stparm(Fr = ilogit(tvFr_logit + nFr_logit)) // Fraction for 1st-order
    stparm(Dur = tvDur * exp(nDur)) // Duration for 0-order
    # ... stparm for Cl, V, Ka ...

    # ... fixef/ranef statements ...
}

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: The sequence is created by using the same parameter (D in this example) to define the duration of the zero-order process and the tlag (time lag) of the first-order process. This elegantly ensures the second process begins precisely when the first one ends. This method is compatible with closed-form (cfMicro) solutions.

Example Code (Sequential ZOFO):

test() {
    # This model uses cfMicro because the time-dependency is handled by
    # the event scheduler (dosepoint), not by changing the ODEs themselves.
    cfMicro(A1, Cl / V, first = (Aa = Ka))
    C = A1 / V

    # SEQUENTIAL PATHWAYS:
    # 1. First-Order Pathway:
    # A fraction 'Fr' of the dose is sent to the absorption depot 'Aa',
    # but its absorption is DELAYED by 'D' time units.
    dosepoint(Aa, tlag = D, bioavail = Fr)

    # 2. Zero-Order Pathway:
    # The remaining fraction '(1-Fr)' is absorbed into the central
    # compartment 'A1' over a DURATION of 'D' time units, starting at time 0.
    dosepoint(A1, bioavail = 1 - Fr, duration = D)

    # Define the fraction 'Fr' to be absorbed first-order (after the lag)
    stparm(Fr = ilogit(tvFr_logit + nFr_logit))
    fixef(tvFr_logit = c(, 0, ))
    ranef(diag(nFr_logit) = c(0.09))

    # Define the duration/tlag parameter 'D'
    stparm(D = tvD * exp(nD))
    fixef(tvD = c(, 1, ))
    ranef(diag(nD) = c(0.09))

    # ... observe/error statements ...
    # ... stparm, fixef, ranef for Cl, V, Ka ...
}

C. Splitting a Single Dose into Different Administration Profiles (e.g., Bolus + Infusion)

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 on the same compartment. The crucial split argument is added to the first dosepoint statement. This allows a single dose amount from one column in the data (e.g., AMT) to be divided between the two statements using their respective bioavail options.

Example (50% Bolus, 50% Infusion over 3 time units):

This model takes a single dose amount and administers 50% of it as an instantaneous bolus into Aa, while simultaneously starting a 3-hour infusion of the other 50% into the same Aa compartment.

PML Code:

test(){
    # Define the central compartment and concentration
    deriv(A1 = - (Cl * C) + (Aa * Ka))
    deriv(Aa = - (Aa * Ka))
    C = A1 / V

    # DOSE SPLITTING:
    # 1. Define the bolus part. 'bioavail=(.5)' takes 50% of the dose.
    #    'split' is ESSENTIAL to allow the same data column to be used for dosepoint2.
    dosepoint(Aa, bioavail = (.5), split)

    # 2. Define the infusion part on the same compartment 'Aa'.
    #    'bioavail=(.5)' takes the other 50% of the dose.
    #    'duration=(3)' defines the zero-order infusion time.
    dosepoint2(Aa, bioavail = (.5), duration = (3))

    # ... observe/error and parameter definition statements ...
    error(CEps = 30)
    observe(CObs = C + CEps)
    stparm(V = tvV * exp(nV))
    stparm(Cl = tvCl * exp(nCl))
    stparm(Ka = tvKa * exp(nKa))
    # ... fixef/ranef statements ...
}

Corresponding Column Definition Mapping:

# Because 'split' is used, both dose() and dose2() can map to "AMT".
dose(Aa<-"AMT")
dose2(Aa<-"AMT")

D. Parallel Absorption with Independent Lag Times

For maximum flexibility, especially with complex formulations, you can assign independent lag times to each parallel absorption pathway. This allows for modeling scenarios where two different release mechanisms from a single dose not only have different profiles (e.g., zero-order vs. first-order) but also start at different times.

Key Concept: Use a separate tlag argument in each dosepoint statement. Each tlag can be linked to a different structural parameter, allowing them to be estimated independently.

Example (Simultaneous ZOFO with Independent Lags): This model splits a dose into a first-order pathway and a zero-order pathway, where each can start after its own unique delay.

PML Code Snippet (Focus on dosepoints):

    # 1. First-Order Pathway with its own lag time 'Tlag_FO'
    dosepoint(Aa, bioavail = Fr_FO, tlag = Tlag_FO)

    # 2. Zero-Order Pathway with its own lag time 'Tlag_ZO'
    dosepoint(A1, bioavail = 1 - Fr_FO, duration = Dur_ZO, tlag = Tlag_ZO)

    # ... parameter definitions for Tlag_FO and Tlag_ZO would be separate ...
    stparm(Tlag_FO = tvTlag_FO * exp(nTlag_FO))
    stparm(Tlag_ZO = tvTlag_ZO * exp(nTlag_ZO))

Related Topics: dosepoint Statement, split, bioavail, duration, tlag

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.

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

End of Part 10

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)

Title: stparm Statement

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

Title: dosepoint Statement

Title: urinecpt Statement

Part 5: Observation and Error Models

Title: observe Statement

Title: error Statement

Title: Error Models

Title: BQL Handling

Title: Finding Extrema (peak function and alternative methods)

Title: Table Output

Part 6: Delays

Title: delay Function

Title: gammaDelay Function

Title: delayInfCpt Statement

Title: transit Statement Keywords: transit compartment, absorption

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

A. Simultaneous (Parallel) Absorption (Zero-Order and First-Order)

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

C. Splitting a Single Dose into Different Administration Profiles (e.g., Bolus + Infusion)

D. Parallel Absorption with Independent Lag Times

Title: Modeling Multiple Elimination Pathways

Part 11: Command-Line Execution

Title: NLME Engine Workflow Overview

Title: Command-Line Syntax and Arguments

Title: Data Input: Primary Dataset (-d1)

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

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

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

Title: Data Input Flags Summary

Title: Command-Line Environment Setup

Key Locations and Variables

How it Works in Practice

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

Estimation Method (-m)

Title: `stparm` 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

Title: `dosepoint` Statement

Title: `urinecpt` Statement

Title: `observe` Statement

Title: `error` Statement

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

Title: `delay` Function

Title: `gammaDelay` Function

Title: `delayInfCpt` Statement

Title: `transit` Statement Keywords: transit compartment, absorption

Title: `if` and Ternary Operator

Title: `sleep` Function

Title: `sequence` Block

Title: `real` and `double` Statements

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...`)