Defining NLME Models

# Defining NLME Models

Nonlinear Mixed Effects (NLME) models are central to pharmacometric modeling. A model is the structure which includes the dynamical equations, the structure of the parameters (the names, domains, and constraints), and the observables. This page documents the two interfaces for defining a NLME model:

1. The @model Domain-Specific Language (DSL) is for simplified definitions of NLME models with standard naming assumptions.
2. The function-based interface is for defining NLME models via Julia functions, allowing for full flexibility and efficiency.

We recommend that all users start with the @model DSL, and computationally-inclined pharmacometricians who are comfortable with more programmatic development (or those who need the enlarged feature-set) may wish to utilize the function-based interface. Both interface with the proceeding simulation and estimation tooling in the same manner.

### Quick Note on Probability Distributions

Many of the NLME model definition portions require the specification of probability distributions. The distributions in Pumas are defined by the Distributions.jl library. All of the Distributions.jl Distribution types are able to be used throughout the Pumas model definitions. Multivariate domains defines values which are vectors while univariate domains define values which are scalars. For the full documentation of the Distribution types, please see the Distributions.jl documentation

## The @model DSL

The @model DSL allows for simplified NLME definitions. This interface can also be used to define simpler linear and probability distribution models without the mixed effects. The components of a model (in order) are as follows:

1. @param defines the fixed effects and other parameters of the model, along with the domains and constraints on the parameters (for estimation).
2. @random (optional) defines the random effects via probability distributions on the parameters.
3. @covariates (optional) defines the covariates of the model.
4. @pre defines the pre-processing collation between the parameters, random effects, and covariates for the definition of the dynamical parameters.
5. @vars (optional) defines aliases that can be used in the proceeding blocks.
6. @init (optional) defines the initial conditions for the dynamical model.
7. @dynamics (optional) defines the dynamical model, either by its differential equation or its analytical solution.
8. @derived defines the derived observables from the dynamical model's solution, including the error distributions.
9. @observed (optional) defines post-processing on the observables sampled from the error distributions.

All of these blocks allow the use of Julia functions defined outside of the macro.

### @param: Parameters

This block defines the structure of the parameters within the model. Parameters are defined by an in (or ∈, written via \in) statement specifying the Domain type that the parameter resides in. For example, to specify θ as a real scalar, one would write:

θ in RealDomain()
θ ∈ RealDomain()

Many of these domains allow specifying bounds, for example, we can specify θ as a scalar between 0.0 and 1.0 via:

θ ∈ RealDomain(lower=0.0, upper=1.0)

For the full specifications of the domain types, please the Domains page. Additionally, parameters can be defined via a probability distributions, in which case the values are defined via ~. For example, we can say that θ comes from a standard normal distribution via:

θ ~ Normal(0,1)

Implicitly, the domain of θ is the support of the probability distributions. Thus for this example, the domain of θ is the same as RealDomain(). However, if the value of θ is not specified in the parameter list during simulation, θ will automatically be sampled from this distribution. Additionally, this probability distribution can be thought of as the prior distribution on θ and is utilized in Bayesian estimation routines.

The parameters block is a list of parameter domain definitions, for example:

@param begin
θ ∈ VectorDomain(3, lower=[0.0,0.0,0.0], upper=[20.0,20.0,20.0])
Ω ∈ PSDDomain(2)
Σ ∈ ConstDomain(0.1)
end

Note that this block is for the structure and not the values of the parameters. The values are defined when invoking simulation or estimation so that they can more easily be modified.

### @random: Random Effects

This block defines the structure of the random effect sampling process. These structures are given by ~ statements to probability distributions which may be defined by parameters. For example, if Ω is a positive-definite matrix, we can specify that η is defined as a sample from a multivariate normal distribution with covariance matrix Ω via the statement:

η ~ MVNormal(Ω)

The @random block is defined by a list of such statements, like:

@random begin
η ~ MVNormal(Ω)
κ ~ MVNormal(Π)
end

### @covariates: Covariates

The @covariates block defines the names of the covariates. This is a simply a list of names, such as:

@covariates wt sex height

or

@covariates begin
wt
sex
height
end

Covariates in the model match the structures they inherit from the data defined in the Subject.

### @pre: Pre-Processing

The @pre block defines the pre-processing collation for the definition of the dynamical parameters from the fixed and random effects. These values are specified by equality (=) statements. For example, one may specify that the parameter Ka is defined by the first value of θ and the first value of η, we can write the command:

Ka = θ[1] * exp(η[1])

Standard Julia syntax can be used within this block, any externally defined Julia functions can be used in this block, and the resulting variables can be any Julia type. One consequence of allowing these values to be any Julia type is that the pre-processed variables can be Julia functions. For example, we can define Ka as a time-dependent function by using Julia's anonymous function syntax:

Ka = t -> t*θ[1]

@pre is defined by a list of such equality statements, for example:

@pre begin
Ka = θ[1]
CL = θ[2]*exp(η[1])
Vc = t -> t*θ[3]*exp(η[2])
end

#### Dosing Control Parameters

Special parameters, such as lag, are used to control the internal event handling (dosing) system. For more information on these parameters, see the Dosing Control Parameters (DCP) page.

### @vars: Variable Aliases

The @vars block defines aliases which can be used in the proceeding blocks. In the @init and @dynamics blocks the statement is interpreted to take place at the current solver time, while in the @derived and @observed the values alias the time series along the solution. An alias is defined by an equality (=) statement. For example, to define conc as an alias for the dynamical variable Central divided by the parameter V, we would write the equation:

@vars begin
conc = Central / V
end

Note that the variable t for time can be utilized within these expressions. Inside the @init and @dynamics blocks it stands for the solver time, while in the @derived and @observed blocks it stands for the current time in the time series. For example,

@vars begin
conc_t = conc / t
end

is the value of conc divided by t.

The @vars block is defined by a list of such equality statements, such as:

@vars begin
conc   = Central / V
conc_t = conc / t
end

Note that the special value := can be used to define intermediate statements that will not be carried outside of the block.

### @init: Initial Conditions

This block defines the initial conditions of the dynamical model in terms of the parameters, random effects, and pre-processed variables. It is defined by a series of equality (=) statements. For example, to set the initial condition of the Response dynamical variable to be the value of the 5th term of the parameter θ, we would use the syntax:

Response = θ[5]

The block is then given by a sequence of such statements, such as:

@init begin
Response1 = θ[5]
Response2 = Kin/Kout
end

where Kin and Kout were defined earlier in the @pre block.

Any variable omitted from this block is given the default initial condition of 0. If the block is omitted, then all dynamical variables are initialized at 0.

Note that the special value := can be used to define intermediate statements that will not be carried outside of the block.

### @dynamics: The Dynamical Model

The @dynamics block defines the nonlinear function from the parameters to the derived variables via a dynamical (differential equation) model. It can currently be specified either by an analytical solution type or via an ordinary differential equation (ODE) (for more types of differential equations, please see the function-based interface).

The analytical solutions are defined in the Dynamical Problem Types page and can be invoked via the name. For example,

@dynamics OneCompartmentModel

defines the dynamical model as the OneCompartmentModel.

For a system of ODEs, the dynamical variables are defined by their derivative expression. A derivative expression is given by a variable's derivative (specified by ') and an equality (=). For example, the following defines the value Depot by it's ODE:

Depot' = -Ka*Depot

where Ka was defined in the @pre block. Variable aliases defined in the @vars are accessible in this block. Additionally, the variable t is reserved for the solver time. For example, if Ka(t) was defined as a function in the @pre block, then the value of Ka at solver time can be utilized in the derivative expression via:

Depot' = -Ka(t)*Depot

This is utilized for handling constructs such as time-varying covariates.

Note that any Julia function defined outside of the @model block can be invoked in the @dynamics block.

### @derived: Derived Observables

The @derived block defines the derived observables of the NLME model. They can be defined by any combination of the parameters, random effects, covariates, preprocessed variables, dynamical variables, and aliases. In this block, the value t is the time series which matches the array given in subject.time. The dynamical variables are an array which matches t in size, where var[i] is the value of the dynamical variable at time t[i]. Any aliases of a dynamical variable are also a time series.

Observables can either be defined by equality statements = or by a distribution with ~. For example, the equality statement

conc = @. Central / V

defines an array conc to be output from the model. Notice that we used Julia's broadcast syntax (@.) for specifying that every value of Central is to be divided by V. Note that any standard Julia syntax (and externally defined functions) are allowed in this block.

Error models are defined by ~ statements to probability distributions. For example, the following defines a time series of Normal distributions centered around the value of conc with a variance dependent on conc and ϵ:

@derived begin
dv ~ @. Normal(conc,conc*ϵ)
end

The likelihood of these distributions are utilized in the maximum likelihood and Bayesian estimation routines. Additionally, values in the @observed block are sampled from these error models.

The @derived block is defined by a list of these expressions, for example:

@derived begin
conc = @. Central / V
dv ~ @. Normal(conc,conc*ϵ)
end

Note that the special value := can be used to define intermediate statements that will not be carried outside of the block.

As a convenience, tie-ins with the included Noncompartmental Analysis (NCA) suite are given with via the @nca macro. For example, we can perform an NCA analysis via:

@derived begin
conc = @. Central / V
dv ~ @. Normal(conc,conc*ϵ)
nca := @nca conc
end

to build an NCASubject using the time series given by the derived or dynamical variable conc. Once defined, the functionality of the Noncompartmental Analysis (NCA) can be used to define derived variables via NCA diagnostics, for example:

@derived begin
conc = @. Central / V
dv ~ @. Normal(conc,conc*ϵ)
nca := @nca conc
auc =  NCA.auc(nca)
thalf =  NCA.thalf(nca)
cmax = NCA.cmax(nca)
end

Notice that the @derived block can mix values of different types (such as arrays and scalars) in the output.

### @observed: Sampled Observations

The @observed block allows one to define output variables based on the sampled values from the error model. These are given by equality statements (=) which can utilize the parameters, random effects, covariates, dynamical variables, and any sampled derived variables. For example, if we had defined:

@derived begin
dv ~ @. Normal(conc,conc*ϵ)
end

then we can add the simulated AUC of the concentration with the error model's stochasticity (the dv values of one simulation) by utilizing the NCA features from within the @observed block as follows:

@observed begin
nca := @nca dv
sampled_auc = NCA.auc(nca)
end

If no @observed block is specified, then the results of a simulation will simply be the derived values and the samples from the error models.

## The PumasModel Function-Based Interface

The PumasModel function-based interface for defining an NLME model is the most expressive mechanism for using Pumas and directly utilizes Julia types and functions. In fact, under the hood the @model DSL works by building an expression for the PumasModel interface! A PumasModel has the constructor:

PumasModel(paramset,random,pre,init,prob,derived,observed=(col,sol,obstimes,samples,subject)->samples)

Notice that the observed function is optional. This section describes the API of the functions which make up the PumasModel type. The structure closely follows that of the @model macro but is more directly Julia syntax.

### The paramset ParamSet

The value paramset is a ParamSet object which takes in a named tuple of Domain types. These Domain types are defined on the Domains page. For example, the following is a value ParamSet construction:

paramset = ParamSet((θ = VectorDomain(4, lower=zeros(4)), # parameters
Ω = PSDDomain(2),
Σ = RealDomain(lower=0.0),
a = ConstDomain(0.2)))

### The random Function

The random(param) function is a function on the parameters. It takes in the values from the param input named tuple and outputs a ParamSet for the random effects. For example:

function random(p)
ParamSet((η=MvNormal(p.Ω),))
end

is a valid random function.

### The pre Function

The pre function takes in the param named tuple, the sampled randeffs named tuple, and the subject data and defines the named tuple of the collated preprocessed dynamical parameters. For example, the following is a valid definition of the pre function:

function pre(param,randeffs,subject)
(Σ  = param.Σ,
Ka = param.θ[1],  # pre
CL = param.θ[2] * ((subject.covariates.wt/70)^0.75) *
(param.θ[4]^subject.covariates.sex) * exp(randeffs.η[1]),
V  = param.θ[3] * exp(randeffs.η[2]))
end

The output can be any valid Julia type. Notice that the covariates are specified via the subject.covariates field.

#### Dosing Control Parameters

Special parameters in the return of the pre function, such as lag, are used to control the internal event handling (dosing) system. For more information on these parameters, see the Dosing Control Parameters (DCP) page.

### The init Function

The init function defines the initial conditions of the dynamical variables from the collated preprocessed values col and the initial time point t0. Note that this follows the DifferentialEquations.jl convention, in that the initial value type defines the type for the state used in the evolution equation.

For example, the following defines the initial condition to be a vector of two zeros:

function init(col,t0)
[0.0,0.0]
end

### The prob DEProblem

The prob is a DEProblem defined by DifferentialEquations.jl. It can be any DEProblem, and the choice of DEProblem specifies the type of dynamical model. For example, if prob is an SDEProblem, then the NLME will be defined via a stochastic differential equation, and if prob is a DDEProblem, then the NLME will be defined via a delay differential equation. For details on defining a DEProblem, please consult the DifferentialEquations.jl documentation.

Note that the timespan, initial condition, and parameters are sentinels that will be overridden in the simulation pipeline. Thus, for example, we can define prob as an ODEProblem omitting these values as follows:

function onecompartment_f(du,u,p,t)
du[1] = -p.Ka*u[1]
du[2] =  p.Ka*u[1] - (p.CL/p.V)*u[2]
end
prob = ODEProblem(onecompartment_f,nothing,nothing,nothing)

Notice that the parameters of the differential equation p is the result value pre.

### The derived Function

The derived function takes in the collated preprocessed values col, the DESolution to the differential equation sol, the obstimes set during the simulation and estimation, and the full subject data. The output can be any Julia type on which map(f,x) is defined (the map is utilized for the subsequent sampling of the error models). For example, the following is a valid derived function which outputs a named tuple:

function derived(col,sol,obstimes,subject)
central = sol(obstimes;idxs=2)
conc = @. central / col.V
dv = @. Normal(conc, conc*col.Σ)
(dv=dv,)
end

Note that probability distributions in the output have a special meaning in maximum likelihood and Bayesian estimation, and are automatically sampled to become observation values during simulation.

### The observed Function

The observed function takes in the collated preprocessed values col, the DESolution sol, the obstimes, the sampled derived values samples, and the full subject data. The output value is the simulation output. It can be any Julia type. For example, the following is a valid observed function:

function observed(col,sol,obstimes,samples,subject)
(obs_cmax = maximum(samples.dv),
T_max = maximum(obstimes),
dv = samples.dv)
end

Note that if the observed function is not given to the PumasModel constructor, the default function (col,sol,obstimes,samples,subject)->samples which passes through the sampled derived values is used.