Tutorial for Modeling Delays in Biological Systems in the NONMEM Software

Abstract

ABSTRACT

Delays in biological systems are a common phenomenon. The models for delays require specialized mathematical and numerical techniques such as transit compartments, delay differential equations (DDEs), and distributed DDEs (DDDEs). Because of mathematical complexity, DDEs and particularly DDDEs are infrequently used for modeling. DDEs are supported by most pharmacometric programs. Recently, DDDEs have been implemented in NONMEM that greatly improve the applicability of this technique in pharmacokinetic and pharmacodynamic (PKPD) modeling. The objective of this tutorial is to provide examples of PKPD models with delays and demonstrate how to implement them in NONMEM. All examples provide a brief description of the biology and pharmacology underlying model equations, explain how they are coded in the NONMEM control stream, and discuss results of data analysis models were used for. NONMEM codes for all models are presented in supporting information (Data S1). The tutorial concludes with a discussion of the pros and cons of presented delay modeling techniques with guidelines for which one might be preferred given the nature of the delay, available data, and the task to be performed.

Introduction

This tutorial demonstrates modeling delays for PK/PD problems in NONMEM. Biological phenomena benefiting from being modeled with delays include: absorption, biophase distribution, indirect response, cell lifespan, or cell maturation. Absorption delay of a drug delivered into a non-IV site is easy to model, and uses the ALAG model parameter available in NONMEM. Delays due to an indirect response can be described with differential equations modeling synthesis and degradation of a receptor or enzyme that then influences a product or cellular activity. An example is the effect of warfarin on prothrombin complex activity (PCA) in blood [1]. Warfarin inhibits vitamin K-dependent synthesis of PCA, which causes a temporal decrease followed by return to the baseline when the drug washes out. Such models offer at least a semi-mechanistic explanation of delay between the time a drug may affect the synthesis, degradation, or activity of the receptor, and the action observed downstream as a biomarker. These target-mediated drug disposition models can be readily modeled using the standard ordinary differential equation tools such as ADVAN6, ADVAN8, ADVAN13, ADVAN14, and ADVAN15. Other delays, such as cell lifespan and cell maturation, can be particularly challenging to model and benefit from modeling methods such as transit compartments or delay differential equations (DDEs) available in NONMEM. For example, the hematopoietic system for red blood cells (RCBs) consists of development from stem cells to intermediate stages to erythrocytes [2]. Production of RCB precursors in bone marrow is stimulated by erythropoietin. The cells remain at each maturation stage for a period before moving on to the next stage to be released to the blood. The average lifespan of a human erythrocyte in blood is about 120 days. Similar maturation delays followed by elimination from circulating blood due to lifespan expiration exhibit platelets simulated by thrombopoietin [3] and neutrophils stimulated by granulocyte-colony stimulating factor [4]. These types of delays will be discussed in this tutorial.

The easiest delay to conceptualize assumes every cell resides at a stage the same amount of time, with little variability. While rarely are biological processes so synchronized, it is reasonable to model in this way if the actual variability of lifespan is no more than, say, about 10% of the average lifespan. Here the simplifying assumption is reasonable, in that residence in that stage is but one of many processes that occur in the biological system that is modeled and may have little impact on the final product measured. Estimating lifespan as a point (delta) distribution is sufficient and allows for a rapid analysis. Such models are most efficiently evaluated by use of DDE solvers, such as ADVAN16, ADVAN17, and ADVAN18, one of the tools covered in this tutorial.

If the lifespan that cells reside in a stage varies by 10% or more, one should consider the average lifespan and its standard deviation. It is usually sufficient to capture the mean and variance of possible lifespans of a stage in a complex biological system, and assume one of two reasonable distributions often used. A typical distribution for biological/cellular processes is the gamma distribution. A typical distribution for lifespans of individual organisms or patients is the Weibull distribution. Along with discrete delay (mean lifespan with zero variance), these two distributions will be used in the examples of this tutorial, although NONMEM can also model a variety of other distributions. A typical process describing a distributed lifespan consists of using the input rate to a state variable (such as a cell stage) as the output rate, but distributed over time. Mathematically, this is done by the method of convolution, where the input stage will be convolved (spread over time) with the gamma or Weibull distribution function, and used as the output rate for the present stage. This process is difficult to model unaided, so various tools in NONMEM are available to facilitate the modeling.

One can model a distributed lifespan process in NONMEM in an approximated manner. The more precise one wishes to model the distribution, the more terms are required, and the more computation time is needed. One method is transit compartments inserted between the input stage and the next stage. A special distribution is the gamma distribution with an integer value n for its shape parameter, called the Erlang distribution, which, if n is not too large, a series of n transit compartments of a common rate constant applied to each transit compartment can be inserted between two stages, to provide the mathematically equivalent (and precise) distribution of integer-n gamma distributed lifespans. If a non-integer n is required, or one wishes to estimate the n, or a non-gamma distribution is required, then a series of transit compartments are produced, and a weighted average of the transit compartments is produced such that the inputted state variable is approximately convolved according to the shape of the desired distribution. This is not trivial, so NONMEM provides a tool that creates the appropriate equations and rate constant expressions that will approximate the desired distribution. This expanded transit compartment method can be solved using standard ODE solvers, such as ADVAN6, 8, 13, 14.

The more accurate the rendering of the distributional delay process required, the more transit compartments (and therefore differential equations) are needed. For the Erlang distribution with integer n, 1/RSTD² transit compartments are needed, where RSTD is the relative standard deviation. So, an RSTD of 10% (0.1) requires 100 simple transit compartments. Hence, once the standard deviation is less than 10%, it is often preferable to approximate the problem with a discrete (point distribution) delay, as described above.

Another method is to convert the convolution integral expression that represents the distributed delay process into an approximate sum of the input at discretely delayed times multiplied by the probability density function of the desired distribution, much like one approximates the integral area under the curve by a weighted sum of the concentrations, as in the trapezoidal rule. Any such approximation of an integral with a weighted sum is called a quadrature process, and NONMEM provides a number of quadrature methods. As with any discrete sum approximating a continuous integral, the greater the number of quadrature points (delay time positions), the more accurate the convolution (integration) process, and of course the more computationally expensive. In addition to the familiar trapezoidal rule, others are available, the most accurate (per quadrature point) being the Gaussian quadrature method, which is the default. Because this method uses discrete time delays, it requires the use of the DDE solvers but does not add differential equations like the general transit compartment method does.

One can obtain exact distributional time-delay representations rapidly calculated relative to the other distributional methods, if the only state variable upon which a time-delay is processed is the final state variable, and does not serve as an input to any other state variable (including its own, such as in a feedback process). Furthermore, the derivative of the state variable is governed by a simple form of input rate minus output rate, in which the output rate is a distributional delayed form of the input rate. One can use the Repetition Variables RPTO/RPTI method that has been available since NONMEM VI. The cumulative distribution function (CDF) is required, and NONMEM has these functions for many distributions. An example will be shown using the Weibull CDF.

The Introduction to NONMEM 7 Guide [5] serves as the reference manual for the implementation of the delay differential methods in NONMEM 7, and this tutorial will refer to it for additional details on available options. In addition, the NONMEM7 Technical guide [6] describes the mathematical motivation for these methods.

Terms specific to time delay modeling and their meaning are listed in Table 1. Table 2 lists various delay methods that will be discussed with their pros and cons, and Table 3 is a decision tree to help you decide what method to choose. Refer to these tables as you are going through the examples.

TABLE 1 Glossary of terms for time delay systems.

Item	Meaning	Placement
TAUy	User defined discrete time delay, defined in $PK, where y = integer. Example: TAU1=EXP(MU_4+ETA(4)) There may be more than one time delay (TAU2, TAU3, etc.)	$PK or $DES before placement of AP_x_y
AP_x_y	Variable defining the past, for state variable A(x), delayed by time TAUy. Defined in $DES. Examples: AP_1_1 = Y0; constant past for A(1), delay TAU1 AP_4_1 = RC0; constant past for A(4), delay TAU1 AP_2_5 = AA*EXP(BB*T); Variable past for A(2), delay TAU5	$DES after TAUy has been defined, but before AD_x_y is used
AD_x_y	State variable A(x) delayed by time TAUy, used in $DES. Example: DADT(1) = KG*(1.0-AD_1_1/YSS)*A(1)	$DES after TAUy and AP_x_y have been defined.
RPTO	Repetition output value, set by user Enables use of repetition feature. Initialize Set repetition value with base 1, for example:$INFNIF (ICALL.EQUATION 0) RPTO=1$PKIF it is not the first record in the subject,IF (NEWIND.EQUATION 2) RPTO=-1This causes $DES to restart integration from time when RPTO was set to 1 (at time 0) to the present TIME	$INFN, $PK
RPTI	Repletion input value, set by NONMEM, for example:IF RPTI=0, then this is not a repeat call Use this to update the delay time in preparation for the repetition, for example:IF (RPTI.EQUATION 0) TI=TIME	$PK
;DDE	Inform ddexpand utility to supply code for a DDE solver, rather than supplying code for a method of steps process	Beginning of control stream, before any NONMEM interpreted code
;NETA	Inform ddexpand the number of etas in the problem, if some of the etas are symbolically defined, such as ETA(CL). Example: ;NETA = 6
;NNy	Number of quadrature positions, or number of transit differential equations (compartments) to be used in convolution of distributed time delay system y (y is integer from 1 to …)	After ;DDE
APT_x_y	Use transit compartment method to provide an Erlang Ersatz Gamma distributed time delay. Variable defining the past, for state variable A(x) (or variable named x if x is not integer), delay system y.	$DES, after TAUy has been defined, but before ADT_x_y is used
ADT_x_y	State variable A(x) or variable x time delayed by distributed delay system y for an Erlang Ersatz Gamma distribution	$DES, after TAUy and APT_x_y have been defined.
;PRCy	Precision for a general distribution delay system y: distribution will be used from time 0 up to tlast = distribCDFINV(1-PRCy). Example: ;PRC1 = 0.03	After ;NNy
;DISTRIBy = name_of_distribution(vector)	Distribution to be used for general transit compartment delay system y, and name of vector to be used for defining parameter values. Example: ;DISTRIB1 = WEIBULLCDF(VG)	After ;PRCy
Define parameters for distribution	Example:VG(2)=NU1VG(3)=TAU1	$PK or $DES
APG_x_y	Provide a general distribution time delay. Variable defining the past, for state variable A(x), delay system y.	$DES, after TAUy has been defined, but before ADG_x_y is used
ADG_x_y	State variable A(x) time delayed by general distributed delay system y	$DES, after TAUy and APG_x_y have been defined.
;DISTRIBy = name_of_distribution(vector, number of parameters)	Distribution to be used for packaged quadrature delay method, vector name to use for parameter inputs to distribution function, and number of parameters to define distribution. Will use discrete time delay calls to DDE solver at NNy time delay positions. Example:;DISTRIB1=WEIBULL(VG,2)	After ;PRCy
;ALGy	Optional. Define quadrature method for specifying the NNy time delay positions and weights of the distributed delay system y. Use Legendre type Gaussian function for quadrature process, which uses a uniform weight process, for example: ;ALG1 = GQ(U)	After ;DISTRIB
;DOCz	Start of quadrature loop z, to allow convolution equations to be built into the control stream, for example: ;DOC1 15 TS CS GQ(G,0.5,TEND) 15 values TS (from 0 to TEND) and associated weights CS, use Laguerre (gamma function of shape 0.5) type Gaussian quadrature	$DES
;ENDDOCz	End of equations that belong in the loop. For example: ;ENDDOC1 ends loop that started with ;DOC1 DOC/ENDDOC loops may be nested. But must have different z indices.	$DES

TABLE 2 Description of various delay methods. Listed are the methods of delay that may be implemented.

Number	Delay method	Description	Pros	Cons
1	Discrete time delay	Models a biological step (such as an intermediate cell type) with a single lifespan value	Easy to model, defining just TAUy, AP_x_y, and AD_x_y. Rapid computation compared to distributed lifespan No additional differential equations needed Fairly rapid calculation	Not realistic to suppose every cell's duration at a given stage is exactly the same: inaccurate if actual lifespan variability is greater than 10% Requires a DDE solver
2	RPTO/RPTI repetition variable method	Uses repetitive reintegration from time 0 to TIME of the data record to model a distributed output rate of a final stage cell type by convolution	Moderate effort to model, requiring monitoring of RPTI and setting of RPTO in a few places, and equation describing the convolution Fairly rapid calculation No additional differential equations needed Convolution is exact to level of TOL, ATOL Any ODE solver will work	Can only be used to model output (removal) of final stage cell type or process. This distributed output rate may not serve as input to another or its own compartment (such as a feedback process)
3	Erlang Ersatz Gamma Distribution	Uses transit compartment to provide a gamma distributed convolution (if NU = integer), or approximate gamma distribution convolution (if NU is not integer)	Easy to model, requiring only TAUy, APT_x_y, and ADT_x_y variables, and ;NNy declaration of maximum number of compartments (maximal NU) For greater than 10% variability of distributed delay times, Fairly rapid calculation Requires addition of fewer than 100 differential equations Can be used for any transitional process, not just end processes Any ODE solver will work	Requires addition of up to 100 differential equations for variability near 10% Can only be used for gamma, or near gamma, distributed time delays FOCE can be problematic because of discontinuity, should use SAEM/IMP/BAYES
4	General distribution transit compartment	Uses a weighted sum of transit compartment state variables to approximate any probability distribution delay	Moderate amount of modeling effort, requiring declaration of number of transit compartments (;NNy), declaration of distribution (;DISTRIBy), in addition to use of APG_x_y, and ADG_x_y in model code If the distribution emulation needs only to be approximate, may need as few as 10 transit compartments Can be used for any distribution and any transitional process, not just end processes Any ODE solver will work	If fidelity to the intended distribution needs to be within 1% or closer, can require 100–200 transit differential equations
5	Packaged Quadrature Delay	Uses a discrete, delay differential equation (DDE) solver to approximate any probability distribution delay	Moderate amount of modeling effort, requiring declaration of number of quadrature points (;NNy), declaration of distribution (;DISTRIBy), in addition to use of APG_x_y, and ADG_x_y in model code Can provide fidelity to intended distribution within 1% for as few as 10–15 quadrature (;NNy) points Can be used for any distribution and any transitional process, not just end processes No differential equations are added to the system Computation time is less expensive than general transit method (4), when fidelity to intended distribution is needed	Requires DDE solver Computation time can be expensive relative to general transit method if fidelity to intended distribution is not required Requires DDE solver Can only convolve state variables A() Or, the Variable's derivative is expressed as a DADT, so its state variable is the rate to be convolved Or, use DAE solver to equate a state variable with the variable
6	Inline Code Quadrature Delay	Uses a discrete delay differential equation (DDE) solver to approximate any second independent variable integration	Can provide fidelity to intended distribution within 1% for as few as 10–15 quadrature (;NN) points Can be used for any distribution and any transitional process, not just end processes Can be used also for nondistribution convolutions, and any variable, not just state variable, may be convolved No differential equations are added to the system Computation time is less expensive than general transit method (4), when fidelity to intended distribution is needed	Considerable amount of modeling effort, requiring in-line coding into control stream file a ;DOC loop construct with equations describing the convolution Requires DDE solver Computation time can be expensive relative to general transit method if high fidelity to intended distribution is not required

TABLE 3 Decision table for selecting method.

Consideration number	Consideration	How to test	If yes	If No
1	Is variability of time delay to the process greater than 10%, or adding variability to the delay process improves the goodness of fit	Use method 3 to model the time delay portion in question and fix the NU value to 10, 30, and 100. Do the resulting OFV's vary by more than 20 units? Or simulate a single subject, set NU to 10, 30, 100, do the predicted values change by a significant amount?	Go on to consideration 2	Use method 1
2	Is the state variable that results from a convolution used to describe the derivative to any other state variable (including its own)	Is the state variable that is a product of the convolution used in any DADT() equation	Go on to consideration 3	Use method 2
3	Is the probability distribution used important	If there is a specific distribution suggested by theory, then the answer is yes Otherwise: Use method 4 using ;NN value of 15–30, then fit the data trying out a Gamma distribution, then by a Weibull distribution, then by lognormal distribution. If difference in OFV is greater than 20 units, answer is yes	Go on to Consideration 4	Use method 3
4	Is a gamma distribution or gamma-like distribution needed, but can be approximated	If there is a specific distribution suggested by theory, and it is not the gamma distribution, then the answer is no Otherwise: If in test for consideration 3 the Gamma distribution provided the lowest OFV, the answer is yes	Use Method 3 if not using FOCE Use method 4 with ;DISTRIB = GAMMA and ;NN < =30 if using FOCE	Go on to consideration 5
5	Is the type of distribution important	If there is a specific distribution suggested by theory, then the answer is yes Otherwise, If the test for consideration 3 resulted in an OFV for one of the tested distributions (you may expand the type of distribution you test) that is significantly lower than the others, then the answer is yes	Go to consideration 6	Use method 3
6	Must the convolution method be true to the intended distribution to within 1% or less	Use method 4 using the specific distribution needed, and try out ;NN = 10, then ;NN = 30, then ;NN = 100. If the OFV improvement from 30 to 100 is important enough to warrant needing ;NN > 30, answer is yes	Go to consideration 7	Use Method 4
7	Can the convolution process be modeled as a distribution convolved with a model variable	If the model requires an integration with a secondary time variable that cannot be represented as a standard convolution, or that is not convolved with a probability distribution, the answer is no	Use method 5 If the model variable to be convolved is not an A() state variable, use $AES, equilibrium compartment method to define a state variable as equivalent to the variable to be convolved	Use method 6

Example 1: Modeling Discrete Delay Times in a Biological Process

We shall consider a simple example of a problem using discrete delay times, which require DDE solvers for greatest efficiency. The pros and cons of the discrete delay method are listed in Table 2 (method 1). For many problems, it is sufficient to model the average time delay of a cellular/biological process using the discrete time delay model, without requiring a distribution with variance, or the data does not have sufficient granularity to allow for a distributed delay model. Or, if the relative standard deviation is less than 10%, which for a gamma distribution would require a NU value of 100 or greater, it is sufficient to model this event with a discrete time delay. Discrete time delay problems, especially when solved with the DDE solvers, do not add as much additional computational burden as distributed delay models, and certainly do not add any additional differential equations. See decision tree Table 3 for the decision process of when it would be suitable to model discrete time delays.

To prepare your control stream for these solvers, select ADVAN16, ADVAN17, or ADVAN18 in the $SUB record.

Consider the logistic growth model that consists of the following logistic equation to model an organism's population growth [7],1$\frac{\mathit{dA}\left(t\right)}{\mathit{dt}}=k_{g}A\left(t\right)\left(1.0-\frac{A\left(t\right)}{\mathit{yss}}\right),A\left(0\right)=A_0$

Only state variable values at the present time t of the integration are needed, and no time delay process is involved. The population growth rate is proportional to the number of members presently $A\left(t\right)$, multiplied by the fraction remaining to reach steady-state yss.

A more realistic rendition of population growth is that each member has a lifespan τ, and the number of individuals dying at the present time t is proportional to the number of individuals that were born a time τ before the present t [8]:2$\frac{\mathit{dA}\left(t\right)}{\mathit{dt}}=k_{g}A\left(t\right)\left(1.0-\frac{A\left(t-\tau \right)}{\mathit{yss}}\right),A\left(0\right)=A_0,-\tau \le t\le 0$

The state variable at time t − τ, $A\left(t-\tau \right)$, requires that the solver store all the values of A() that were generated in the past, from time 0 up to the present time t, and use those values at will, in accordance with model requirements. Because t − τ can be negative at times t < τ, a definition of A() is required for times before the numerical integration process started at t = 0. Usually this is a constant past that is equivalent to some steady-state condition before numerical integration began at time t, such as A0, the starting state variable value at time 0. This model, along with several other discrete delay models, have been implemented in NONMEM and described in [9].

The above logistic problem is an easy, one-derivative problem that will initiate us into the method of DDEs in a population analysis setting.

For the NMTRAN control stream, we define certain reserved variables that represent A(t − τ), so NONMEM can insert the appropriate state variable value at time t − τ.

Three types of reserved variables are defined for a DDE problem, where τ is represented by TAU (also listed in Table 1):

TAUy:

User defined time delay, defined in $PK, where y = integer. For example: TAU1=EXP(MU_4+ETA(4))

There may be more than one time delay (TAU2, TAU3, etc.)

AP_x_y:

Variable defining the past, for state variable x, delayed by time TAUy. Defined in $DES.

Examples:

AP_1_1 = Y0; constant past for A(1), delay TAU1

AP_4_1 = RC0; constant past for A(4), delay TAU1

AP_2_5 = AA*EXP(BB*T); Variable past for A(2), delay TAU5.

AD_x_y:

State variable A(x) delayed by time TAUy, used in $DES.

Example:

DADT(1) = KG*(1.0-AD_1_1/YSS)*A(1)

The rate of change of A(1) depends on the A(1) at the present time T, and on AD_1_1, which is the A(1) at T-TAU1.

Consider the base code without time delay components:$DESDADT(1)=KG*(1.0-A(1)/YSS)*A(1)

The user modifies the code by inserting definitions of TAU1 and utilizing AP_1_1, and replacing A(1) with AD_1_1 as appropriate:

Code	Explanation
$PK
TAU1=EXP(MU_4+ETA(4))	Define TAU1, with random effect as needed
$DES
AP_1_1=Y0	;Define Past for A(1), in this case steady-state pre-dose level
DADT(1)=KG*(1.0-AD_1_1/YSS)*A(1)	Use AD_1_1 in equations as state A(1) delayed by time TAU1

We will use RADAR5 DDE SOLVER ADVAN16, which is the most efficient DDE solver.

EM methods (importance sampling, SAEM) are the most efficient estimations for DDE:

$EST METHOD=IMP INTERACTION MCETA = 10

MCETA adds a random search that makes EM (as well as FOCE) analysis more robust. Importance sampling creates random samples of parameters (via etas), and weights their significance according to goodness of fit to the individual's data and population fixed effects. A weighted average of these random samples of parameters is used to update the population parameter (theta).

MU referencing maps which random samples of etas are to be used to update which thetas.

EM methods are efficient when the to-be-estimated thetas are MU-referenced.

MU_x is defined as the typical value of a population parameter. It represents the arithmetic mean of the normal distribution of individual parameters associated with the population parameter:

TVCL=THETA(2)	Typical value of clearance
LTVCL=LOG(THETA(2))	Log of typical value of clearance
MU_2=LTVCL	MU_2 is the log of the typical value of clearance
CL=EXP(MU_2+ETA(2))	MU_2 must be associated arithmetically with ETA(2)

The full Example 1 control stream is listed in Code list 1 (example1\ho1.ctl, bolded where different from usual: Code List items are in the Code_List.docx document with the Data S1 as a convenience to have all control stream code collated in one document). When executing the control stream, all on one line, use the -dde option:

..\nmfe76 ho1.ctl ho1.res -prdefault –dde

This control stream will evaluate the problem using the importance sampling method.

The -dde option enlists the ddexpand utility program to modify the control stream and add lines to define the AD_x_y variables and make special calls to the RADAR5 delay algorithm. The ddexpand utility sets up arguments to and calls a built-in function DDEFUNC, which will return the value of state variable A() from time T-TAU1 earlier and use this value where AD_x_y was positioned by the user. You can see the final control stream submitted to NMTRAN in ho1.ctl_dde if you are curious. Results to this analysis are shown in Figure 1 and Table 4.

[IMAGE OMITTED. SEE PDF]

TABLE 4 Population analysis results for delayed logistic growth model, with a discrete time delay.

Parameter	KG	Y0	YSS	TAU1
Mean	0.199 ± 0.003	0.976 ± 0.016	9.89 ± 0.153	4.82 ± 0.081
Omega STD	0.0929 ± 0.0123	0.0877 ± 0.0120	0.0842 ± 0.0108	0.0821 ± 0.0128
Residual STD	0.0568 ± 0.0008
-2LL	−889.6

Example 2: Modeling Discrete Delay Times Between Stages in a Biological Process 2

A more complicated example of discrete delay times consists of erythropoietin action on RCB maturation, with time delays (lifespan) for each stage [10] (Figure 2A). This model describes the pharmacokinetics of recombinant human erythropoietin and its effect on RBCs and hemoglobin (Hb) in rats [11]. An analysis of the effect of an erythropoietin-mimetic product on RCB production in healthy volunteers using NONMEM's discrete delay system was published earlier [9].

[IMAGE OMITTED. SEE PDF]

The time course of erythropoietin (EPO) concentration in plasma (C) following an injection IV bolus D_IV is quantified by a target-mediated disposition model:3$\frac{\mathit{dC}}{\mathit{dt}}=k_{\mathrm{EPO}}-k_{\mathrm{on}}·R·C+k_{\mathrm{off}}·\mathrm{RC}-\left(k_{\mathrm{el}}+k_{\mathrm{pt}}\right)·C+\frac{k_{\mathrm{tp}}·A_T}{V_p}$4$\frac{d{A}_T}{\mathit{dt}}=k_{\mathrm{pt}}·C·V_p-k_{\mathrm{tp}}·A_T$5$\frac{\mathit{dR}}{\mathit{dt}}=k_{\mathrm{syn}}-k_{\mathrm{on}}·R·C+k_{\mathrm{off}}·\mathrm{RC}-k_{\mathrm{deg}}·R$6$\frac{\mathit{dRC}}{\mathit{dt}}=k_{\mathrm{on}}·R·C-\left(k_{\mathit{off}}+k_{\mathit{int}}\right)·RC$

Endogenous EPO is synthesized at constant rate k_EPO into compartment C, and exogenous drug is administered into the same compartment intravenously. EPO distributes to peripheral compartment (A_T) at first-order rates k_pt and k_tp and is eliminated from plasma at first-order rate k_el. EPO binds to and dissociates from EPO receptors (R) at rates k_on and k_off, respectively. Free EPO receptors are synthesized at zero-order rate k_syn and degraded at first-order rate k_deg. Upon binding to the receptor EPO forms drug-receptor complex RC that is internalized at first-order rate k_int. EPO stimulates the production of RBC precursor cells in bone marrow P₁ and P₂. The precursors P₁ maturate to precursors P₂ after the time T_P1. The precursor cells P₂ are released to the circulation as reticulocytes RET after time T_P2. Reticulocytes become mature RCBs RBC_M after time TRET, and RBC_M are removed from the circulation do to senescence after their lifespan T_RBC. The rate of the removal of cells from the pools P₁, P₂, RET, and RBC_M is determined by their lifespans. The resulting system of DDEs describing the states RET and RBC_M is as follows:7$\begin{array}{cc}\frac{d\mathrm{RET}}{\mathit{dt}}=\hfill & k_{\text{in}}·S\left(t-T_{P1}-T_{P2}\right)·S\left(t-T_{P2}\right)·I\left(t-T_{P1}-T_{P2}\right)-\hfill \\ \hfill & k_{\text{in}}·S\left(t-T_{P1}-T_{P2}-T_{\mathrm{RET}}\right)·S\left(t-T_{P2}-T_{\mathrm{RET}}\right)·I\left(t-T_{P1}-T_{P2}-T_{\mathrm{RET}}\right)\hfill \end{array}$8$\begin{array}{cc}\frac{d\mathrm{RB}{\mathrm{C}}_{\mathrm{M}}}{\mathit{dt}}=\hfill & k_{\text{in}}·S\left(t-T_{P1}-T_{P2}-T_{\mathrm{RET}}\right)·S\left(t-T_{P2}-T_{\mathrm{RET}}\right)·I\left(t-T_{P1}-T_{P2}-T_{\mathrm{RET}}\right)-\hfill \\ \hfill & k_{\text{in}}·S\left(t-T_{P1}-T_{P2}-T_{\mathrm{RET}}-T_{\mathrm{RBC}}\right)·S\left(t-T_{P2}-T_{\mathrm{RET}}-T_{\mathrm{RBC}}\right)·I\left(t-T_{P1}-T_{P2}-T_{\mathrm{RET}}-T_{\mathrm{RBC}}\right)\hfill \end{array}$

Since states P₁ and P₂ do not appear in Equations (7) and (8) there is no need to include them in the final model aiming at describing the time courses of RET, RBC, and Hb. RBCs are the sum of RET and RBC_M:9$\mathrm{RBC}\left(t\right)=\mathrm{RET}\left(t\right)+\mathrm{RB}{\mathrm{C}}_{\mathrm{M}}\left(t\right)$and Hb is derived from RBCs by means of the mean corpuscular Hb concentration MCH.10$\mathrm{Hb}\left(t\right)=\mathrm{MCH}·\mathrm{RBC}\left(t\right)$

The zero-order production rate for the precursor cells P₁ is multiplied by the stimulatory function S(t) defined by the drug–receptor complex plasma concentration RC:11$S\left(t\right)=1+\frac{S_{\max }·\mathrm{RC}\left(t\right)}{{\mathrm{SC}}_{50}+\mathrm{RC}\left(t\right)}$where S_max is the maximal stimulation and SC₅₀ is the concentration of RC eliciting 50% of the maximal stimulation. The same stimulatory function multiplies the production rate for P₂. The model accounts for a negative feedback of an increase in Hb$\mathrm{\Delta Hb}\left(t\right)=\mathrm{Hb}\left(t\right)-\mathrm{Hb}\left(0\right)$on the production of P1. The inhibitory function I(t)$I\left(t\right)=1-\frac{I_{\max }·\mathrm{\Delta Hb}\left(t\right)}{{\mathrm{IC}}_{50}+\mathrm{\Delta Hb}\left(t\right)}$multiplies the production rate constant k_in. Here, I_max is the maximal inhibition and IC₅₀ is the increase in the Hb concentration eliciting 50% of the maximal inhibition. The system is at steady state prior to EPO dose injection at time t = 0:12$C\left(t\right)=C_0,\text{for}t<0,\text{and}C\left(0\right)=C_0+\frac{D_{\mathrm{IV}}}{V_p}$13$A_T\left(t\right)=\frac{k_{\mathrm{pt}}·C_0·V_p}{k_{\mathrm{tp}}},\text{for}t\le 0$14$R\left(t\right)=R_0,\text{for}t\le 0$15$\mathrm{RC}\left(t\right)=R{C}_0,\text{for}t\le 0$16$\mathrm{RET}\left(t\right)=\mathrm{RE}{\mathrm{T}}_0,\text{for}t\le 0$17$\mathrm{RB}{\mathrm{C}}_{\mathrm{M}}\left(t\right)=\mathrm{RB}{\mathrm{C}}_0-\mathrm{RE}{\mathrm{T}}_0,\text{for}t\le 0$where C₀, R₀, RC₀, and RET₀ are the baseline (steady-state) values for C, R, RC, and RET, respectively. RBC₀ is the baseline value for RBC. To ensure that the system steady state is at baseline values, the following relationship between the remaining model parameters needs to be imposed:18$\mathrm{R}{\mathrm{C}}_0=\frac{k_{\mathrm{on}}·R_0·C_0}{k_{\mathrm{off}}+k_{\mathrm{int}}}$19$k_{\mathrm{EPO}}=k_{\mathrm{el}}·C_0+k_{\mathrm{int}}·R{C}_0$20$k_{\mathrm{syn}}=k_{\mathrm{deg}}·R_0+k_{\mathrm{int}}·R{C}_0$21$\mathrm{RE}{\mathrm{T}}_0=\frac{T_{\mathrm{RET}}·\mathrm{RB}{\mathrm{C}}_0}{T_{\mathrm{RET}}+T_{\mathrm{RBC}}}$22$k_{\text{in}}=\frac{\mathrm{RE}{\mathrm{T}}_0}{T_{\mathrm{RET}}·{\left(1+\frac{{\mathrm{S}}_{\max }·{\mathrm{RC}}_0}{{\mathrm{SC}}_{50}+{\mathrm{RC}}_0}\right)}^2}$

Model parameters used for simulation are listed in Table 5. The variables pertaining to delayed processes are as follows.

TABLE 5 Parameter values obtained from [10] and used for the erythropoietin stimulated red blood cell production model of Example 2.

Parameter	Value
DIV, pmol/kg	92.51
Vp, mL/kg	56.94
kel, h−1	0.2256
kpt, h−1	0.2092
ktp, h−1	0.1721
kint, h−1	0.8228
kdeg, h−1	0.1133
kon, pM−1 h−1	0.01132
koff, h−1	1.297
R0, pM	63.2
C0, pM	3.248
RBC0, 1012 cells/L	6.128
MCH, 10−1 pg/cell	2.0
TP1, h	42.97
TP2, h	3.02
TRET, h	72.33
TRBC, h	1440
S max	3.48
SC50, pM	1.7
I max	1.0
IC50, pM	1.79

TP1=THETA(17)	Time to mature from precursor 1 to 2
TP2=THETA(18)	Time to mature from precursor 2 to reticulocytes
TRET=THETA(19)	Time to mature from reticulocytes to circulating RBC
TRBC=THETA(20)	Lifespan of RBC in circulation

Six composite lifespans are needed for the model, using the above base life spansTAU1=TP1+TP2 TAU2=TP2 TAU3=TP1+TP2+TRET TAU4=TP2+TRET TAU5=TP1+TP2+TRET+TRBC TAU6=TP2+TRET+TRBC

The ddexpand program uses the AP_x_y as a placeholder to insert code that initializes arguments and calls DDEFUNC that returns the state variable value of A(x) from a time TAUy before the present integration time t (at t-TAUy) and places this value in the variable AD_x_y of each distinct y. Here are some examples of the AP_x_y's defined for this model:

AP_4_1=RC0	Past for State A(4), using TAU1 as delay: RC(T-TP1-TP2)
AP_4_2=RC0	Past for State A(4), using TAU2 as delay: RC(T-TP2)
AP_5_5=RET0	Past for State A(5), using TAU1 as delay: RET(T-TP1-TP2-TRET-TRBC)
AP_6_5=RBCM0	Past for State A(6), using TAU5 as delay: RBCM(T-TP1-TP2-TRET-TRBC)

Some variables set up for easy expression of differential equations, and use of delay states: AD_x_y:X1 = 1+SMAX*AD_4_1/(SC50+AD_4_1)&#x2026;I5 = 1&#x2212;IMAX*(MCH*(AD_5_5+AD_6_5)&#x2212;HB0)/(IC50+(MCH*(AD_5_5+AD_6_5)&#x2212;HB0))

Differential Equations using variables that were defined with delay states&#x2026;DADT(5)=KIN*X1*X2*I1&#x2212;KIN*X3*X4*I3DADT(6)=KIN*X3*X4*I3&#x2212;KIN*X5*X6*I5

The complete control stream is shown in code list 2 (example2\epodd.ctl). Execute the problem as:

nmfe76 epodd.ctl epodd.res -prdefault –dde

and the problem will run. See epodd.ctl_dde in your run directory for the processed control stream after the execution.

Parameters used are listed in Table 5, and RBC profile is in Figure 2B. To assist ddexpand, add the following commented lines, preferably at the beginning of the control stream (see Table 1):

;DDE	Indicate to ddexpand that system is to be expanded for use by DDE solver
;NETAS=6	Number of etas in problem: Required by ddexpand if some etas are symbolically defined (such as ETA(CL))

These techniques will also be handy when we use the discrete differential equation solvers to solve distributed delay problems, shown in subsequent sections.

Example 3: Solving Distributed Delayed Output Rate Problems Using Repetition Variables RPTO/RPTI (nmVI)

Conditions may occur requiring a delay modeled in a distributed manner, particularly for delay times distributed by more than 10% relative to the average time delay (see decision tree Table 3, for method 2). In one class of problems, the derivative state variable is governed by a simple form of input rate minus output rate, in which the output rate is a delayed form of the input rate. The distribution of the delay times is governed by a probability density function l, so this can be mathematically described as23$\frac{\mathit{dN}\left(t\right)}{\mathit{dt}}=k_{\text{in}}\left(t\right)-{\int }_0^{\infty }l\left(s\right)k_{\text{in}}\left(t-s\right)\mathrm{ds}$

The second term on the right hand side is the convolution integral, and it weights input rate delayed by time s, k_in(t − s), with pdf function l(s), and sums over all possible time delays, from s = 0 to s = infinity. This differs from the discrete delay which we considered earlier, where a single delay time s = tau is considered. This distributed delay problem can be solved using the RPTO/RPTI repetition variable method (Method 2 listed in Table 2). The output is delayed by a t − s, with possible values of s varying from 0 to infinity, distributed according to l(s). We simplify with further derivation, resulting in:24$\frac{\mathit{dN}\left(t\right)}{\mathit{dt}}=\left(1-L\left(\mathrm{TI}-t\right)\right)k_{\text{in}}\left(t\right)$where L() is the CDF of l(), and TI is the TIME given in the data record, so that when the ODE solver has completed the integration up to time TI, the result is25${\int }_0^{\mathit{TI}}\mathit{dN}\left(t\right)={\int }_0^{\mathit{TI}}\left(1-L\left(\mathrm{TI}-t\right)\right)k_{\text{in}}\left(t\right)\mathit{dt}=N\left(\mathrm{TI}\right)$

NONMEM outputs $N\left(\mathrm{TI}\right)$ as the predicted value for the value of TIME of that data record. The detailed derivation is provided in [6].

We shall use the cumulative WEIBULL distribution for L for this example. The Weibull density is as follows:

For $\alpha \in R^{+}$, $\beta \in R^{+}$, $y\in R^{+}$,26$\text{WEIBULL}\left(y|\alpha ,\beta \right)=\frac{\alpha }{\beta }{\left(\frac{y}{\beta }\right)}^{\alpha -1}\exp \left(-{\left(\frac{y}{\beta }\right)}^{\alpha }\right)$where $\alpha$ is the shape parameter, and $\beta$ is the scale parameter, in units of time. The mean transit or delay time for a Weibull distribution is27$\mathrm{MTT}=\mathit{\beta \Gamma }\left(1+1/\alpha \right)$

For $\alpha \ge 1$, average time is approximately the scale parameter $\tau$. The relative standard deviation (STD/Mean Time) is:28$\text{RSTD}=\sqrt{\frac{\Gamma \left(1+2/\alpha \right)}{\Gamma {\left(1+1/\alpha \right)}^2}-1}$

As shape parameter $\alpha$ increases in value, the RSTD decreases, and distribution is concentrated more and more narrowly around the mean time. The CDF is the integral of the density:29$\begin{array}{c}L\left(\text{integral}\text{from}0\text{to}\mathrm{y}\right)=\text{WEIBULLCDF}\left(y|\alpha ,\beta \right)=1.0-\exp \left(-{\left(\frac{y}{\beta }\right)}^{\alpha }\right)\hfill \end{array}$

And30$\begin{array}{c}1-L\left(\text{integral}\text{from}0\text{to}\mathrm{y}\right)=\exp \left(-{\left(\frac{y}{\beta }\right)}^{\alpha }\right)\hfill \end{array}$

We convolve the Weibull distribution on an input function governed by an Smax relation to the drug level in the central compartment (example3\smax_weibull_rep_foce.ctl), linked with a Michaelis–Menten PK problem from [12], shown in Code List 3:

c = A(1)/V	Concentration in central compartment
kin = kin0*(1.0+(Smax*c)/(SC50+c))	Input rate
if(ti-­t>=0.0)DD=EXP(−((TI-­T)/BB)**AA)	DD=1-L(t)=1-WEIBULLCDF(TI-T), AA=shape parameter, BB=scale parameter
DADT(1) = -(Vmax*A(1))/(V50+A(1))
DADT(2) = DD*(KIN-KIN0)	=dN(t)/dt=(1-L(TI-t))*(kin-kin0)

For examples in this tutorial that use the Weibull distribution, the shape parameter will be alternately labeled as AA or NU, and the scale parameter will be alternately labeled as BB or TAU.

Note the RPTO part of the code is in $INFN and $PK.$INFNIF (ICALL.EQ.0) RPTO=1 ;enables use of repetition feature$PKIF (RPTI.EQ.0) TI=TIMEIF (NEWIND.EQ.2) RPTO=-1

Because TI = TIME changes with each data record, NONMEM needs to repeat the integration from t = 0 to t = TI using the TI = TIME value appropriate for the present record. The RPTO/RPTI setup causes repeated integration from T = 0 to T = TIME = TI of the present record, and for each record, thus performing a convolution of KIN-KIN0, with the integrated result A(2) = N(t) being the total number of cells at time t, a result of an output rate that is equal to the production rate Weibull distribution delayed.

Figure 3A shows the shape of the Weibull PDF for NU = 7.0, Tau = 6.5, and Figure 3B shows the effect of the convolution of the kin with the Weibull distribution.

[IMAGE OMITTED. SEE PDF]

The RPTO convolution process can be implemented only if the derivative of the state variable as a function of time can be represented in a k-k*l format, where k is the input/production rate. This convolution method using the repetition variables RPTO/RPTI will not work for general convolution problems. Suppose we have, as in the first problem:DADT(1)=KG*(1.0-(A(1)*l(t))/YSS)*A(1)

This expression cannot be put into a k − k*l format, and therefore cannot be treated with a (1 − L)*k type of solution, which is the only type that the RPTO/RPTI method can solve properly. For general distributed problems, we need the distributed delay methods described later, which provide a second variable of integration process, by creating an extended partitioning via transit compartments, or by discrete time partitions using the DDE solver.

The NONMEM repetition variables method produces a distributional convolution of kin exactly calculated up to TOL and ATOL, whereas the distributed transit or distributed quadrature methods described in the next sections only approximate these convolutions. Furthermore, any ODE solver can be used with this method, and a specialized DDE solver is not required. The repetition variables method does add a small amount of computation time, in that the numerical integration by the ODE solver is repeated from T = 0 to T = TIME of each data record. The more data records in the data set, the greater the computation time. However, as related above, the RPTO/RPTI method cannot be used for general convolution problems and cannot be used for more general distributed delay problems, whereas distributed transit and distributed quadrature delay methods, shown in the next section, can be used to solve more general convolution problems.

Example 4: Transit Compartments for Gamma Distributed Time Delays (Erlang Ersatz Gamma Distributed Delay Method)

The gamma distribution can be emulated by the simple transit compartment method. This transit compartment method is listed in Table 2 as method 3, along with its pros and cons. See Table 3 to determine under what conditions this method would be suitable for your needs, such as > 10% variability of time delay, and when an approximate gamma distribution is appropriate. The gamma density is:

For $\nu \in R^{+}$, $k\in R^{+}$, $y\in R^{+}$,31$\mathit{\text{GAMMA}}\left(y|\nu ,k\right)=\frac{k^{\nu }}{\Gamma \left(\nu \right)}{y}^{\nu -1}\exp \left(-\mathit{ky}\right)$where $\nu \in R^{+}$ is the shape parameter, and $k\in R^{+}$ is the scale parameter in units of 1/time. The mean transit time is:32$\mathrm{MTT}=\nu /k$

The relative standard deviation (STD/mean time) is:33$\text{RSTD}=\frac{1}{\sqrt{\nu }}$

As the shape parameter $\nu$ increases in value, the RSTD decreases, and the distribution is concentrated more and more narrowly around the MTT.

Emulation of the gamma distribution is done by adding several ODE's to the NMTRAN control stream, so that state variable A(x) is filtered through $\nu$ (NU) transit compartments in which the rate constant of transfer KTR (k) is in one direction, much as is done when filtering a dose input through a series of transit compartments. The KTR_y is related to the mean transit time (MTT) of TAUy and NUy transit compartments as:

k = KTR_y = NUy/TAUy

This is modeled in the same way as discrete delay, but instead of specifying AP_x_y and AD_x_y, you specify APT_x_y and ADT_x_y (T for transit). The transit compartment method is equivalent to a gamma distribution when NUy is an integer (called the Erlang distribution). If NUy is a non-integer, the ddexpand utility will evaluate CEILING(NUy) transit compartments and provide the linear weighted average of the CEILING(NUy) transit compartment and the (CEILING(NUy)-1)th transit compartment. For example, if NUy = 2.75, then 0.25 of state variable to transit compartment 2 and 0.75 of transit compartment 3 will be added together and serve as state variable ADT_x_y. While this is strictly speaking not the gamma distribution of NUy = 2.75, the difference is typically only a few percent in the time shift (i.e., time at which the same output value is reached for the true gamma versus this approximation is a few percent), and is very easy and fast to calculate. We shall call this the Erlang Ersatz Gamma distribution. Consider the following user-created nmtran control stream (example4\smax_gammat_imp.ctl), listed fully in Code List 4:

;DDE	Inform ddexpand that code is to be expanded for use with a DDE solver
;NN1=10	Maximum number of transit compartments (NN1 > =NU1)
$PK
…
TAU1=EXP(MU_3+ETA(3))	Declaration of average time delay (MTT) 1
NU1=EXP(MU_4+ETA(4))	Shape parameter of Gamma distribution
$DES
c = A(1)/V
kin = kin0*(1.0 
+(Smax*c)/(SC50+c))
APT_KIN_1=KIN0	Past for variable to be convolved, KIN, with time delay system 1
DADT(1) = 
−(Vmax*A(1))/(V50+A(1))
DADT(2) = kin - 
ADT_KIN_1	KIN is delayed by time system 1, using GAMMA(x, NU1, NU1/TAU1)

NN1 specifies the maximum number of compartments for expansion, and TAU1 defines the average time delay. The APT_x_y must be defined for transit compartments, even though the transit compartment system does not use past information (other than constant past), because ddexpand uses the APT_x_y entry to identify the type of delay as the simple transit compartment system. Notice also that for the APT/ADT system, despite the letter A beginning its name (ADT_), the variable that is being convolved does not have to be strictly an indexed state variable A(x), but can refer to a variable named x, if x is not an integer.

The call of ddexpand results in the expansion in the DES section shown in Code List 4b. In the $PK block, the following is assessed.

NU_1=NU1	NU1 is NU for delay system 1
XNU_1=MAX(INT(NU_1)+1,0)	XNU_1 = nearest upper integer to NU1
YNU_1=NU_1-XNU_1+1.0	Fractional part of NU1

Suppose the NU1 in the present iteration is 4.67. Then, XNU_1 = 5, and YNU_1 = 0.67.

An integer IX_1 increments along with the computation of each transit compartment, until IX_1 = XNU_1 = 5, and then the most recent two state variables are summed, in this case:IF(IX_==XNU_1) THENSUM_1=T_O*((1.0-YNU_1)*A(6)+YNU_1*A(7)+APT_KIN_1)ENDIF

You can see that SUM_1 used as the final “output compartment” to the transit system is the weighted arithmetic average (using YNU_1 and (1-YNU_1) as the weights) of the two compartments encompassing NU1, the gamma shape parameter.

This method is straightforward in that it uses the transit compartment technique that is familiar to many, and a specialized DDE solver is not necessary to implement. The ddexpand utility merely adds in the transit compartment equations based on the instructions from; NN1 = 10, the APT_1_1, and ADT_1_1. Furthermore, as NU1 is treated as a continuous variable, it can be included in estimation. The resulting parameters compared to the values that were simulated are shown in Table 6, along with the very accurate RPTO method. One caveat, the Erlang Ersatz Gamma distribution is continuous with respect to NU, but its first derivative is discontinuous with respect to NU, and an FOCE analysis can be problematic, as it is with this example.

TABLE 6 Nonlinear mixed effects results for the Smax/SC50 example 4.

	Smax	SC50	TAU	NU	KIN0 (fixed)	OMEGA RSV	SIGMA RSV
Simulated with true Gamma function (RPTO method) ..\example4\smax_gamma_rep_sim.ctl	0.5	0.1	4.5	2.5	190	0.09	0.01
Estimated with true Gamma function (RPTO method) ..\example4\smax_gamma_rep_imp.ctl	0.478 ± 0.019	0.111 ± 007	4.65 ± 0.19	2.59 ± 0.14	190	0.053 to 0.13	0.0101 ± 0.0004
Estimated with Erlang Ersatz Gamma ..\example4\smax_gammat_rep_imp.ctl	0.479 ± 0.019	0.119 ± 0.007	4.67 ± 0.19	2.69 ± 0.14	190	0.049 to 0.12	0.0101 ± 0.0004

Figure 4A shows the difference in the pdf shape between the true Gamma pdf and the Erlang Ersatz gamma pdf. Figure 4B shows the cellular state A(2) when its output function is convolved with the true Gamma or the Erlang Ersatz Gamma distribution. You can see that results are similar, differing by a few percent.

[IMAGE OMITTED. SEE PDF]

It may be noted here that a gamma distribution with NU1 = 1 is an exponential function that can be emulated with a single transit compartment, and if the distribution acts on the state variable itself, it is equivalent to one transit compartment removing the amount in that state by the rate constant ktr (the scale parameter to the exponential distribution). For example, c = A(1)/V APT_2_1=0.0 kin = kin0*(1.0+(Smax*c)/(SC50+c)) DADT(1) = -(Vmax*A(1))/(V50+A(1)) DADT(2) = kin - kout*ADT_2_1

is equivalent to: c = A(1)/V NU1=1 ktr=NU1/TAU1 kin = kin0*(1.0+(Smax*c)/(SC50+c)) DADT(1) = -(Vmax*A(1))/(V50+A(1)) DADT(2) = kin - kout*A(3) DADT(3) = ktr*(A(2)-A(3))

Example 5, 6: General Transit Compartment Method: Creating Transit Compartments for Time Delays of Any Distribution

A general method for distributed delays in which any distribution function may be used is the general transit compartment method, developed by Koch and Schropp [13]. This is listed as method 4 in Table 2. If emulation of the desired distribution does not require high precision (10%–20% imprecision is acceptable) to the intended distribution, then this could be a suitable method that will not add too many differential equations to the problem and execute reasonably quickly (see Table 3 for helping you decide if method 4 is an appropriate approach).

The method calculates a weighted sum of transit compartment states, with weights determined according to the distribution it is to emulate. The derivation of how these weights are calculated is provided in [13].

Instead of using APT_x_y for the past, and ADT_x_y for the delayed state, the user should define in $DES the variables APG_x_y and ADG_x_y, respectively (G for general distribution), for mean transit time TAUy and state variable A(x), and the ddexpand utility will make the expansions.

Consider the logistic example for a Weibull distribution. The user-written control stream may be like that of Code List 5 (example5\logdt_pop4.ctl). Notice the beginning comment lines:

;DDE	Let ddexpand know that a dde solver is to be used
;NN1=10	Number of transit compartments to use. The higher the NN1, the more accurate
;PRC1=0.03	Precision: distribution will be used from time 0 up to tlast=WEIBULLCDFINV(1-PRC1)
;DISTRIB1=WEIBULLCDF(VG)	distribution type, and argument array name to use

The number of transit equations to expand is NNy, for time delay y. Next, the precision PRCy is specified. If not specified, then 1/NNy precision is assumed. If ;PRC is specified without the y index, then the one specified by the most recent NNy is used. Finally, the distribution function is specified (;DISTRIBy, again, with y optional, and by default that specified by the most recent NNy), along with an arbitrary vector name to use with it; in this case, the vector name is VG. The user may use any of the following distribCDF functions, as defined in [5, 6]:

WEIBULLCDF

GAMMACDF

LOGNORMALCDF

GOMPMAKECDF (Gompertz-Makeham)

WEIBULLCCDF (combined Weibull)

UNIFORMCDF

An approximation of the convolution of the distribution with the input state variable A(x) is to be calculated. The integral of the CDF is to be evaluated at evenly spaced time intervals from s = 0 to s = TTEND = distribCDFINV(1-PRC1), spaced apart by delta-s = TTEND/NN1. Thus, imprecision is about (MAX(1/NN1, PRC1)). The lower the PRC1, and the higher the number of transit compartments (NN1) in the expansion, the more precisely the convolution process is calculated. To minimize computational expense and have a rough approximation to the distribution, one could select NN1 = 10 and PRC1 = 0.01.

The user defines appropriate parameters to the distribution, into the vector name allocated to that distribution. In the above example, the vector is VG. As these CDF functions are structured for use with $ABBR FUNCTION (such as described for the probability densities listed in [5]), indices starting at position 2 should have the appropriate parameter values specified (the first element VG(1) is for the delay time variable of convolution integration, and is filled in by ddexpand, as needed). The parameter locations for the distributions are ordered according to the descriptions of these distributions in [5].

$DES records are similar to the APT/ADT and AP/AD in the earlier examples:

$DES
VG(2)=NU1	Shape argument to Weibull
VG(3)=TAU1	Scale argument to Weibull
APG_1_1=Y0	Define past value (when t < 0)
DADT(1)=KG*(1.0-ADG_1_1/YSS)*A(1)	Differential equation

While we can identify the TAU1 parameter as approximately the MTT for NU1 > 1 for the Weibull distribution, the exact MTT is:$\mathrm{MTT}==\mathrm{TAU}1\Gamma \left(1+1/\mathrm{NU}1\right)$as noted earlier.

When the above control stream is submitted to ddexpand (using the -dde option to the nmfe script), we get the resulting control stream listed in Code List 5b: (example5\logdt_pop4.ctl_dde), which adds additional differential equations that effectively carry out the afore-mentioned convolution of distribution with the intended state variable, with the resolution of NN1, and up to the time determined by PRC1, in the $DES record.

Notice that in the ddexpand-ed code NN1 = 15 transit compartments A(2) through A(16) have a common rate constant KK_1 just as in a standard transit compartment process, but they are inputted into the routine TRANSITG_SDIM1 via the VG() array, which returns a weighted sum of these compartments in accordance with the WEIBULL CDF WEIBULLCDF(). In this way, one obtains an approximation of the Weibull distributed convolution of the source compartment A(1).

In theory, one can use any of the CDF versions of the distributions listed in [5], but the CDF inverse functions (distribCDFINV) to many of these have not yet been constructed, which is required for this process to assess the best TTEND to serve as the upper limit of the convolution integral, based on the precision PRC1 requested.

Sometimes one must impose a distributed delay on a general variable, rather than a state variable A(*). To do this, set x = variablename for APG_x_y and ADG_x_y. For example, kin is to be convolved with a Weibull distribution, as in example6\smax_weibulldt_imp4.ctl (code list 6):;DDE;NN1=200;PRC=0.015;DISTRIB=WEIBULLCDF(VG)$DESVG(2)=AA ; shape parameterVG(3)=BB ; scale parameterc = A(1)/Vkin = kin0*(1.0+(Smax*c)/(SC50+c))APG_KIN_1=KIN0DADT(1) = -(Vmax*A(1))/(V50+A(1))DADT(2) = kin &#x2013; ADG_KIN_1&#x2026;

The ADG_KIN_1 is the result of convolution of kin with pdf of Weibull. You can optionally use a probability density function:;DDE;NN1=200;PRC=0.015;DISTRIB=WEIBULL(VG)

but generally results are better when using the CDF.

Figure 4C compares the Output state variable A(2) when using the general transit compartment method and compares with the “true” RPTO method.

Example 7, 8, 9: The Quadrature Delay Method: Using Discrete Delay Solvers for Time Delays of any Distribution, Using the Ddexpand Program and Built-In DDEFUNCG Routine

So far, transit compartments were used to create distributed delays. One can use discrete DDE algorithms to perform a distributed delay as well, which create discrete delay evaluations of earlier times on state variables A(), instead of using the transit compartment method. The description of this method is listed in Table 2 as method 5, along with its pros and cons. The advantage of using the DDE solvers for distributional delays is that no additional differential equations are added to the user's model system, and particularly with built-in DDE solvers ADVAN16–ADVAN18, they are more efficient than the transit compartment distributed method. Furthermore, this method can produce convolution to the intended distribution with higher precision with less computational expense than the general transit compartment method 4, if you require imprecision less than 10% (see Table 3 to decide when method 5 is suitable). The user models the discrete delay solver method in the same manner as distributed delays via transit compartments shown in example 6, with only one change required, and that is to specify the number of parameters to the distribution in the comment line of ;DISTRIB, for example:;DISTRIB=WEIBULLCDF(VG,2)

You must use one of the solvers that are capable of discrete delay solving (ADVAN13, ADVAN14, ADVAN15, ADVAN16, ADVAN17, ADVAN18). Doing so will allow the use of discrete DDE solvers on NN1 at selected times (determined by the quadrature method), weighted according to the CDF of the distribution, using by default a Gaussian quadrature method. The principles of such quadrature techniques for integration are described in [6]. Alternatively, state variables will be assessed on discrete NN1 quantiles spaced from 0 to 1 if ;ALG = GQ(U) is requested. The advantage to this method is that there will not be NN1 differential equations added to the control stream. Instead, a built-in routine called DDEFUNCG will perform the calculations for the NN1 discrete delay states, and no added memory (except for ADVAN18) or added differential equations are needed. For ADVAN18, you may need to add $SIZES MAXNRDS = x, where x is the total delay times, or NN1 in this example.

Let's consider example5 above, and we change just the lines from;DDE;NN1=15;PRC1=0.03;DISTRIB1=WEIBULLCDF(VG)

to

;DDE	Let ddexpand know that a dde solver is to be used
;NN1=10	Number of transit compartments to use. The higher the NN1, the more accurate
;DISTRIB1=WEIBULLCDF(VG,2)	Distribution type, and argument array name to use, plus the number of parameters distribution
;ALG=GQ(U)	Use Legendre type Gaussian function for quadrature process, which uses a uniform weight process

In this example, we request a uniform Gauss quadrature on 10 quantile positions of the WEIBULL distribution.

This is now example7\logdd_pop2.ctl. The ddexpand utility will process the control stream (focusing on the $DES record) that is listed in Code List 7: example7\logdd_pop2.ctl_dde.

Instead of creating additional differential equations in the control stream, the work is completed by the built-in routine DDEFUNCG(). This routine will call the DDE solver at NNx discrete times in the past and weight the state variable values at these times based on the quadrature method used. By default, the ddexpand routine will establish the most suitable quadrature method for performing the convolution integration. Other quadrature methods may be chosen, listed in [5], such as the uniform (Legendre) quadrature method, but generally the default method is appropriate: the Laguerre Gaussian quadrature, known to be the most efficient method for many functions in terms of obtaining greatest precision for the given number of nodes (NNx) requested.

One down-side to using the built-in DDEFUNCG routine and non-DEA discrete delay solvers is variables like KIN cannot be directly delayed and require access via an equivalent state variable A(). For example, you would add something like:

DADT(x) = {expression for derivative of kin with respect to time}.

In the example above, kin could be coded by the user as an additional state variable, A(3), as in Code List 8: example8\smax_weibulldd_imp3.ctl.

We can alternatively use the equilibrium compartment method (using a DAE dde solver), but now apply it for distributed delay (Code list 9, example9\smax_weibulldae_imp.ctl).

Note that the $AES record allows an equilibrium equation to be defined:E(3)=A(3)-kin

The E(3) is calculated by the DAE solver to be always 0, so that effectively.$\mathrm{Kin}=\mathrm{A}\left(3\right)$and the state variable A(3) is always equal to kin, and ADG_3_1 is then kin delayed according to delay time system 1 (which is defined as Weibull distribution delayed).

The biological output Figure 4D shows that Gaussian quadrature on discrete delays is very similar to the true convolution process performed by the RPTO method, differing by less than 0.1%. The computation is severalfold faster than that of the general transit compartment method at NN = 200, and is more accurate.

In this example, the ALG = GQ(U) algorithm was most efficient because it used Gaussian quadrature positions between 0 and 1, and then established time positions according to the quantile values of the Weibull pdf. For this example, it was effective because, as you see in the figure on the Weibull pdf for NU = 6.8 and TAU = 6.5, there is an initial period of time where the pdf is very small and does not become dominant until about time = 3.0. So quadrature positions were strategically located in the regions where the pdf is most dominant. This does not take into account, however, the dominant positions suitable for the state variable that it is convolving upon, so this may or may not always be more efficient than the default positioning by Gaussian quadrature of time values along the pdf. The PRC is not used for ALG = GQ(U), but it is used in the default ALG (when ALG is not specified). Often placing PRC = 1.0D-10 allows coverage of nearly the entire distribution region and is the most suitable when using Gaussian quadrature on time positions. An NN of 10–50 usually provides a precision of at least 1%. The best way to determine these values is to create a single subject simulation using the approximate population parameters, set NN1 to 150, and create table outputs of the predictive function values. Use these as “true” values, then by trial and error, adjust the options (particularly NN) and find the lowest NN that is not more than 1% of the “true” NN = 150 values. Then perform the population analysis using that NN (using NN = 150 often would make the estimation time too long).

Example 10: General Quadrature Integration Using Ddexpand, and Routines GQ and NC

In the previous section, we used a built-in routine called DDEFUNCG that efficiently evaluates quadrature integration of distributed time delays using prepackaged distribution functions. The advantage to this method is that the DDEFUNCG routine retains up to 1000 distribution values that needed to be calculated once for every delta-tau position, and the values were then stored for repeated use. This allows for rapid evaluation. Furthermore, the coding of the convolution process is done for you. The down-side is that (1) only convolution integrations using prepackaged distributions could be used as the kernel function, (2) the past (APG_x_y) could only be constant, and (3) the function being convolved needed to be defined in some way as a state variable A(x) (such as creating a state variable that represents kin in the previous section by specifying dA(x)/dt = dkin/dt, or using equilibrium compartment equations and solving a DAE, such as E(x) = kin). A more versatile method that can perform any quadrature integration of any equation set and has none of these limitations is for the user to code within NMTRAN the specific equations to be integrated, using the DOWHILE feature, along with calling the quadrature routines (GQ or NC) directly (method 6, Table 2).

Revisiting the example8\smax_weibulldd_imp3.ctl problem, we noticed that the state variable A(3) needed to be created that ultimately represented kin, since DDEFUNCG can only utilize time-delay state variables convolved with distribution functions. We can alternatively code the necessary integration more directly within NMTRAN, as shown in this example, example10\smax_weibull_doc_imp.ctl (Code List 10), but requiring the modeler to code more of the equations:

$DES
…
;DOC 15 TS CS GQ(G,0.5,TEND)	Start of quadrature loop, returning 15 values TS (from 0 to TEND), and y values (weights) CS(), use Laguerre (gamma function of shape 0.5) type Gaussian quadrature
TAU1=TS	TAU1 = TS is the delay, and the 15 TS values are incrementally presented one at a time for each cycle of the loop
AP_1_1=0.0	Code placement for discrete delay call, and defining past value
LL=AA/BB*((TS/BB)**(AA-1.0))*EXP(−((TS/BB)**AA))	WEIBULL(TS, AA, BB) TS = delay time variable, AA = shape parameter, BB = scale parameter
CDELAY=AD_1_1/V	State variable A(1) delayed by time TS, divided by V, =C(T-TS)
FF=KIN0*(1.0+(SMAX*CDELAY)/(SC50+CDELAY))	KIN, input rate value
SUM=SUM+CS*LL*FF	Summing convolution element L(TS)*C(T-TS), weighted by Gaussian quadrature weight CS
SUMW=SUMW+CS*LL	Summing only the distribution L(TS), weighted with Gaussian weight CS
;ENDDOC	End of equations that belong in the loop
ADC=SUM/SUMW	With loop 0f 15 cycles done, SUM contains full integral of the convolution L*KIN by Gaussian quadrature. SUMW is full integral of the Weibull density by Gaussian quadrature, and is very nearly 1, but not quite, so this denominator is a minor numerical corrector.
DADT(2)=KINW-ADC	Use the final result ADC as the output rate, which is really the input rate convolved, or smeared, in a Weibull distributed manner over a range of time

The ;DOC and: ENDDOC comment lines are interpreted by ddexpand (remember to have the –dde option switch on the nmfe7x script command line) as follows:

;DOCy tells ddexpand that this is the beginning of DOWHILE loop y (if y is not present, it is assumed 1), and it ends at ;ENDDOCy. The items to be listed on the ;DOCy line are (in order):

• Number of points for the integration, which is 15 in this example.
• Variable of quadrature integration, which is TS in this example.
• Coefficient of quadrature integration, which is CS in this example.
• Quadrature method to be used, along with arguments. GQ or NC may be selected, with arguments to the quadrature method to be supplied as given in intro7.pdf.

Notice how the variable of integration TS and the coefficient of integration CS are used in accumulator variables SUM and SUMW. The FF variable is kin calculated using the state variable A(1) delayed by TS (notice how TAU1 needs to be defined and associated with TS), LL is the Weibull density evaluated at time TS. The FF*LL is therefore the convolution summand, weighted with CS, the quadrature integration coefficient. Although we are using Gauss-Laguerre, with gamma function weighting, the gamma function weighting is removed from the CS, so that the user can replace it with whatever density or kernel function LL is needed (in this case the Weibull density).

It is not essential to normalize the result by dividing by SUMW, as SUMW should be nearly 1 if LL is properly a pdf, but it sometimes helps in accuracy.

The integral need not be a classic convolution with distribution, and more complicated convolutions may be applied. The Lympho example described in [1] shows how the lymphocyte trafficking model may be applied. Also, Krzyzanski and Bauer [14] (Appendices 6 and 9, Data S1) show how ;DOC method may be applied to assessing mean transit times in age-structured cell populations, and show the considerable versatility of this method.

Discussion: A Summary of Considerations for DDEs

Data S1 provide results of the population analysis of the problems presented here. Table 3 lists the decision tree for deciding on the delay method, and how one may test for various conditions in determining the most suitable method. The discrete delay method can be very efficient, and when DDE solvers are involved, requires no additional differential equations, as with examples 1 and 2. If a distributed delay model improves the goodness of fit significantly, or the delay process has a large variance (greater than 10% RSTD) in distribution of possible time delays, but the distribution is unimportant, then the most rapid assessment for distributed delays can be done by using simple transit compartments of fewer than 100, as in example 3.

If the distributional shape of probable time delays is important, then choose a classic distribution such as gamma, Weibull, or lognormal, defined by two parameters. If the only states influenced by distributional time delays are terminal ones (no other states use it as an input, including itself), then the RPTO method can be used with high precision and computational speed, such as example 3. If, however, there are states that are modeled with distributional delays which influence their own or other states, then use the general and efficient discrete quadrature system, such as examples 7–10. Often 10–15 (NN) past positions are required to obtain accurate output with imprecision of 1% or lower. The RTOL/ATOL to the DDE solver should be set to 10 or greater. If you do not have available a DDE solver, then you can use the general transit compartment method that provides a weighted average of the transit compartments such that the results represent approximately the desired distribution, such as example 5 or 6. You can experiment with the number of additional equations required by creating simulations, increasing the number of transit compartment ODE's until the change in the output is at an acceptably low level.

Distributed delay problems often benefit from parallel computing because of the added computation expense. Table 7 provides some of the resources required for the various delay problems. Table 8 compares the estimation results of the various algorithms on the Smax/Weibull model for importance sampling analysis. Note that the NU and its omegas for the transit method are skewed relative to the results of RPTO and quadrature delay.

TABLE 7 Performance diagnostics for various delay methods. Below is a list of the performance times for the various discrete and distributed delay methods described in this paper.

Example	Delay technique	Estimation method	Additional ODE's needed	Time needed per subject (s) per iteration or function call	NN
Example3\logte.ctl	Linear interpolation Erlang (proportional transit – ADT_, erlang ersatz gamma)	Single subject	5	0.0257	5
Example1\ho1.ctl	Discrete delay (AD_)	ITS	0	0.127	0
Example4\smax_gammat_its.ctl	Erlang Ersatz Gamma	ITS	5	0.138	5
Example4\smax_gammat_imp.ctl	Erlang Ersatz Gamma	IMP	5	0.163	5
Example3\smax_weibull_rep_its.ctl	Weibull/RPTO	ITS	0	0.234	0
Example3\smax_weibull_rep_imp.ctl	Weibull/RPTO	IMP	0	0.246	0
Example7\logdd_pop2.ctl	Weibull/PGQD (packaged Gaussian quadrature delay, using function DDEFUNCG)	ITS	0	0.289	10
Example5\logdt_pop4.ctl	Weibull/transit (ADG_)	ITS	15	0.35	15
Example8\smax_weibulldd_imp3.ctl	Weibull/PGQD	IMP	0	2.31	15
Example9\smax_weibulldae_imp.ctl	Weibull/PGQD using DAE	IMP	0	2.45	15
Example10\smax_weibull_doc_imp.ctl	Weibull/IGQD (inline coded Gaussian quadrature delay)	IMP	0	2.54	15
Example3\smax_weibull_rep_foce.ctl	Weibull/RPTO	FOCE	0	4.23	0
Example6\smax_weibulldt_imp4.ctl	Weibull/transit compartment (ADG_)	IMP	200	41.6	200

TABLE 8 Comparison of population parameter estimates for various methods.

Parameter	Smax	Sc50	NU(AA)	Tau(BB)	OM(smax)	OM(sc50)	OM(NU)	OM(tau)	Sigma
Simulated	0.5	0.1	7.0	6.5	0.09	0.09	0.09	0.09	0.001
RPTO (ex3)	0.479 ± 0.019	0.109 ± 0.0085	6.84 ± 0.52	6.61 ± 0.34	0.0752 ± 0.0152	0.0880 ± 0.0508	0.0966 ± 0.0346	0.125 ± 0.026	0.00979 ± 00040
Gen. Transit (ex6)	0.479 ± 0.019	0.107 ± 0.0081	8.31 ± 0.84	6.53 ± 0.34	0.0756 ± 0.0153	0.0939 ± 0.0469	0.165 ± 0.0608	0.126 ± 0.026	0.00980 ± 0.0040
Quad. delay (ex9)	0.479 ± 0.019	0.108 ± 0.0084	6.84 ± 0.53	6.61 ± 0.34	0.0753 ± 0.0152	0.0932 ± 0.0532	0.104 ± 0.0385	0.125 ± 0.026	0.00981 ± 0.00040

In general, use the EM methods over the FOCE method for greater speed.

Conflicts of Interest

R.J.B. is a paid employee of ICON plc., and is the present developer of NONMEM. The other author declares no conflicts of interest.

References

W. J. Jusko and H. C. Ko, “Physiologic Indirect Response Models Characterize Diverse Typyes of Pharmacodynamic Effect,” Clinical Pharmacology and Therapeutics 56, no. 4 (1994): 406–419, https://doi.org/10.1038/clpt.1994.155.

Y. Yan and W. Krzyzanski, “Quantitative Assessment of Minimal Effective Concentration of Erythropoiesis‐Stimulating Agents,” CPT: Pharmacometrics and Systems Pharmacology 2 (2013): e62, https://doi.org/10.1038/psp.2013.39.

W. Krzyzanski, P. Wiczling, P. Lowe, et al., “Population Modeling of Filgrastim PK‐PD in Healthy Adults Following Intravenous and Subcutaneous Administrations,” Journal of Clinical Pharmacology 50 (2010): 101S–112S, https://doi.org/10.1177/0091270010376966.

Y.‐M. Wang, W. Krzyzanski, S. Doshi, J. J. Xiao, J. J. Perez Ruixo, and A. T. Chow, “Pharmacodynamics‐Mediated Drug Disposition (PDMDD) and Precursor Pool Lifespan Model for Single Dose of Romiplostim in Healthy Subjects,” AAPS Journal 12, no. 4 (2010): 729–740, https://doi.org/10.1208/s12248‐010‐9234‐9.

S. L. Beal, L. B. Sheiner, A. J. Boeckmann, and R. J. Bauer, NONMEM 7.6.0 Users Guides (ICON plc, 1989–2024), https://nonmem.iconplc.com/nonmem760.

R. Bauer, Technical Guide on Various Methods in NONMEM 7. NONMEM 7.6.0 Users Guides (ICON plc, 1989–2024), https://nonmem.iconplc.com/nonmem760.

N. Bacaër, “Verhulst and the Logistic Equation,” in A Short History of Mathematical Population Dynamics (Springer, 2011), https://doi.org/10.1007/978‐0‐85729‐115‐8_6.

G. E. Hutchinson, “Circular Cause Systems in Ecology,” Annals of the New York Academy of Sciences 50 (1948): 221–246.

X. Yan, R. Bauer, G. Koch, J. Schropp, J. J. P. Ruixo, and W. Krzyzanski, “Delay Differential Equations Based Models in NONMEM,” Journal of Pharmacokinetics and Pharmacodynamics 48, no. 6 (2021): 763–802, https://doi.org/10.1007/s10928‐021‐09770‐z.

R. J. Bauer, G. Mo, and W. Krzyzanski, “Solving Delay Differential Equations in S‐ADAPT by Method of Steps,” Computer Methods and Programs in Biomedicine 111 (2013): 715–734, https://doi.org/10.1016/j.cmpb.2013.05.026.

S. Woo, W. Krzyzanski, and W. J. Jusko, “Target‐Mediated Pharmacokinetic and Pharmacodynamic Model of Recombinant Human Erythropoietin (rHuEPO),” Journal of Pharmacokinetics and Pharmacodynamics 34 (2007): 849–868, https://doi.org/10.1007/s10928‐007‐9074‐0.

H. Y. Cheng and W. J. Jusko, “Mean Residence Time Concepts for Pharmacokinetic Systems With Nonlinear Drug Elimination Described by the Michaelis‐Menten Equation,” Pharmaceutical Research 5, no. 3 (1988): 156–164, https://doi.org/10.1023/a:1015960806202.

G. Koch and J. Schropp, “Distributed Transit Compartments for Arbitrary Lifespan Distributions in Aging Populations,” Journal of Theoretical Biology 380 (2015): 550–558, https://doi.org/10.1016/j.jtbi.2015.06.018.

W. Krzyzanski and R. J. Bauer, “Pharmacodynamic Age Structured Population Model for Cell Trafficking,” Journal of Pharmaceutical Sciences 113 (2024): 257–267, https://doi.org/10.1016/j.xphs.2023.10.040.