101 Utility Programs
1 table_quant and table_resample Utility Programs for Analyzing $COV SIR Data (NM74)
The sub-section Importance Sampling of the Variance-Covariance of the Parameter Estimates (NM74) in section I.58 $COV: Additional Options and Behavior describes how importance sampling may be used to obtain samples around the minimum of the FOCE/Laplace objective function. After random samples with WEIGHT values is recorded in the .ext file, these results may be further analyzed using two utility programs.
1.1 table_quant
The utility ..\util\table_quant will transform the results in the raw output file, utilizing the WEIGHT column, into a table file with frequencies and cumulative values. The utility is to be executed on the command line as follows:
table_quant root.ext root.qnt delimiter start endwhere delimiter is that used in the input file root.ext (s for space, default), and start and end are the range of iterations to be quantized (default is all non-negative iterations). Note that, even Bayes results, which do not have the WEIGHT column, can be processed with this utility, but the weight is then assumed constant among all samples.
The resulting file, root.qnt, will contain for each item, the sorted value, its frequency (freq), and its quantile position or cumulative probability (cum). One can readily plot the quantile value against quantile position for a cdf plot, or sum(freq) by quantile value range for frequency plot. The R script quantplot.R, or Splus script, quantplot.ssc, available in Pdx-Pop 5.2, can be used to view histograms and cdf plots from these result files, and produce quantile tables, linearly-interpolating the results at the most interesting quantile positions (0.025, 0.5, 0.975, etc). To use these plotting scripts, make sure you modify the header information to point to the desired file name and extension (by default these scripts look for a qnt extension). Also, make sure a .qnt file contains just one table of information.
1.2 table_resample
An alternative to viewing the $COV/SIR results via its weight information is to perform a resampling of the information in the .ext file, with samples weighted according to the WEIGHT column (if no WEIGHT column is present, WEIGHT is assumed to be equal among all samples). This is done using the utility table_resample:
table_resample root.ext root_new.ext delimiter newsize SEED start endwhere delimiter is that used in the input file root.ext (s for space, default), and start and end are the range of iterations to be quantized (default is all non-negative iterations). In addition:
- seed=0: non-randomized expansion of the samples, based on
WEIGHTcolumn - seed>0: randomized starting at seed, with repeated samples allowed
- seed<0: randomized starting at
abs(seed), with repeated samples not allowed
If the user chooses seed=0, then newsize samples will be generated, each line of the original root.ext file being repeated in proportion to its WEIGHT value, and these repeated samples will be placed in root_new.ext. Thus the weight of each sample is physically expressed in the manner of repeated rows of that sample. To assure that integer truncation does not render the smaller weighted samples to be not at all expressed, newsize should be something like 10000, or even 100000. The resulting file, root_new.ext, will have the same structure as a BAYES result file, without the WEIGHT column, and the R script bayesplot.R or Splus script bayesplot.ssc, available in Pdx-Pop 5.2, may be used to view histograms, quantile plots, and quantile tables.
If the user chooses seed>0, then newsize samples will be generated randomly and with replacement, in proportion to the WEIGHT column, and placed in the root_new.ext file. If
newsize>min(oldsize of original file, end-start+1)then it would make sense to choose this option. The seed<0 should be used only if
newsize<<min(oldsize of original file, end-start+1)that is, you just want to pick a few samples.
2 table_compare Utility Program (NM72)
The utility program table_compare will compare the numerical values between two table files produced by the NONMEM $TABLE record, and the user may specify the tolerance for the comparison. The syntax is:
table_compare mytable1.tab mytable2.tab , myprecision.xtl >mydifferences.txtwhere delimiter is {, t s} for {comma tab space}, and myprecision.xtl is a precision specification or control file. Default delimiter is space and default control file is table_compare.xtl.
table_compare mytable1.tab mytable2.tab , S myprecision.xtl >mydifferences.txtIn the above example, the first file is comma delimited, and the second one is space (S) delimited.
If a second character is given to a delimiter, then this is for detecting a continuation marker at the end of a line that is to be continued. If a third character is given as a delimiter, this for detecting a continuation marker at the beginning of the continuing line. Some examples are:
table_compare mytable1.tab mytable2.tab ",&" "S&" myprecision.xtl >mydifferences.txt(double quotes may be needed for DOS commands). In the above example, the first file is delimited by commas between column items, and an & at the end of a line breaks the record across multiple lines. The second file is delimited by spaces between column items, and an & breaks a record across multiple lines.
table_compare mytable1.tab mytable2.tab ",&c" "S&c" myprecision.xtl >mydifferences.txtIn the above example, the first file is delimited by commas between column items, and an & at the end of a line breaks the record, with a c at the beginning of the next line. The second file is delimited by spaces between column items, and an & at the end of a continuing line, and a c at the beginning of the next line.
table_compare mytable1.tab mytable2.tab ",&" "SSc" myprecision.xtl >mydifferences.txtIn the above example, the first file is delimited by commas between column items, and an & at the end of a line breaks the record. The second file is delimited by spaces between column items, and no special character at the end of a continuing line (the S serves as a place-holder for line continuation markers, since apace is too ambiguous as a continuator) and a c at the beginning of the next line.
It is useful to redirect difference results to a file, in this example mydifferences.txt. For example, the user may desire that only relative differences greater than 0.01 be reported. A very simple control file could be:
$PRECISION
ALL=0.01,0.003stating that all columns be compared with a relative difference of 0.01, and absolute difference of 0.003. Precision criteria for specific columns in the tables may also be given:
$PRECISION
ALL=0.01,0.003 WRES=0.1,0.2
CL=0.05,0.02The equation for comparison is, if
ABS(X-Y)>R*MAX(ABS(X),ABS(Y))+Athen the difference is reported, where R is relative difference tolerance, and A is absolute difference tolerance.
3 table_to_xml Utility Program (NM72)
The utility table_to_xml program in the NONMEM ..\util directory can be used to convert additional NONMEM output tables produced during the $EST step into XML formatted files. The syntax is as follows, as an example:
table_to_xml my_results.cov my_results_cov.xml ,where the delimiter may be ,, t, or s for comma, tab, or space. Default delimiter is space. The rules (schema, document type definition) by which the xml file is constructed are given in tables.xsd and tables.dtd, which are in the ..\run or ..\util directory.
table_to_xml my_results.cov my_results_cov.xml ",&c"specifies that the table file may have line continuator characters & and c, as described in the table_compare section.
4 xml_compare Utility Program and its Use for Installation Qualification (NM72)
The utility program xml_compare will compare the contents of two NONMEM report XML files that are produced by NONMEM. The syntax to the command line is:
xml_compare myresult1.xml myresult2.xml myprecision.xtl >mydifferences.txtwhere myprecision.xtl is a precision specification or control file. Default delimiter is space and default control file is xml_compare.xtl. It is useful to redirect difference results to a file, in this example mydifferences.txt.
The control file can be quite elaborate, but it allows specification of various precision values for the many different types of values in the NONMEM report XML file, and to ignore certain entries as well. An example xml_compare.xtl file is in the ..\util directory, and has the following contents:
$IGNORE
monitor
elapsed_time
datetime
covariance_status
termination_status
nonmem(version)
parallel_est
parallel_cov
$PRECISION
GENERAL=0.2,0.2 OBJ_BAYES=2.0,0.0 OBJ_SAEM=0,100.0 OBJ_ITS=0,5.0
EIGENVALUES=2.0,0
OBJ_IMP=0,10.0 OBJ_F=0,5.0
DIAG=0.3,0 OFFDIAG=0,0.5 COR=0.0,0.3 VAR=0.3,0.1 COV=-1.0
OBJ_DIRECT=0,100.0
correlation_o=-1.0 INVCOVARIANCE_O=-1 INVCOVARIANCE_D=-1
etashrinksd=0,20 epsshrinksd=0,10 ebvshrinksd=0,20
etashrinkvr=0,20 epsshrinkvr=0,10 ebvshrinkvr=0,20
METHOD=DIRECT ALL=-1
METHOD=SAEM epsshrinksd=0,20The $IGNORE record will ignore all elements with the substrings that are listed, or just a specific attribute of an element, such as nonmem(version).
Under the $PRECISION record, a
GENERAL=R,Acan be given for most items, where R is the relative tolerance, and A is the absolute tolerance. Following the GENERAL specification, tolerances may be specified for other items.
Two items of identical element and attributes are compared between the two files, where the equation for comparison is, between value X of xml file 1 and value Y of xml file 2,
ABS(X-Y)>R*MAX(ABS(X),ABS(Y))+AThe OBJ_BAYES is given a special test, as it has a standard deviation with it:
STD(X,Y)=SQRT(STD(X)2+ STD(Y)2)
ABS(X-Y)>R*STD(X,Y)+AIn the above example OBJ_BAYES=(2,0) means that if the Bayes objective functions in the two files differ by more than 2 standard deviations, then the difference is noted. Please note that while the above test is suitable for tolerance comparison in an installation qualification setting, this is not an appropriate statistical test for model comparisons.
To ignore an item for comparison, specify -1. To specify an exact comparison, use 0,0. To refer to a particular optimization method, then enter METHOD=SAEM for example, and thereafter, all entries of items pertain to that estimation method, until METHOD is changed. The METHOD attribute may have one of the following settings:
FOCE, ITS, IMP, SAEM, DIRECT, BAYES (for standard and NUTS)4.1 Precision item reference
The total list of items, and their scope, are as follows (R/2 = 1/2 of relative error):
| NAME | DESCRIPTION | DEFAULT (R,A) |
|---|---|---|
| GENERAL | Default to most non-matrix items | 0.2,0.2 |
| DIAG | Diagonal elements of OMEGA/SIGMA estimates |
0.1,0 |
| OFFDIAG | Off-diagonal elements of OMEGA/SIGMA estimates |
0.0,0.2 |
| VAR | Diagonals of variance of estimates | 0.2,0 |
| COV | Off-diagonals of covariance of estimates | 0,0.2 |
| COR | Correlations | 0,0.2 |
| TABLE | Table items listed in NONMEM report file. | GENERAL |
| OBJ_BAYES | BAYES objective function |
1,0 |
| OBJ_SAEM | SAEM objective function |
0,100 |
| OBJ_ITS | ITS objective function |
0,2 |
| OBJ_IMP | IMP/IMPMAP objective function |
0,5 |
| OBJ_DIRECT | Direct sampling objective function | 0,100 |
| OBJ_F | FO/FOCE/Laplace objective function |
0,0.5 |
| EIGENVALUES | Eigenvalues | 2,2 |
| ETABAR | Etabar | GENERAL |
| ETABARSE | Etabar Se | GENERAL |
| ETABARPVAL | Etabar Pval | GENERAL |
| ETASHRINKSD | Eta shrinkage, SD type | GENERAL |
| EPSSHRINKSD | EPS shrinkage, SD type | GENERAL |
| EBVSHRINKSD | ETA Empirical Bayes Variance shrinkage, SD type | GENERAL |
| ETASHRINKVR | Eta shrinkage, variance type | GENERAL |
| EPSSHRINKVR | EPS shrinkage, variance type | GENERAL |
| EBVSHRINKVR | ETA Empirical Bayes Variance shrinkage, variance type | GENERAL |
| RELATIVEINF | ETA Relative information | GENERAL |
THETA |
THETAs | GENERAL |
| OMEGA_D | Omega diagonals | DIAG |
| OMEGA_O | Omega off-diagonals | OFFDIAG |
| SIGMA_D | Sigma diagonals | DIAG |
| SIGMA_O | Sigma off-diagonals | OFFDIAG |
| OMEGAC_D | Omega correlation diagonals | DIAG (R/2,A) |
| OMEGAC_O | Omega correlation off-diagonals | COR |
| SIGMAC_D | Sigma correlation diagonals | DIAG (R/2,A) |
| SIGMAC_O | Sigma correlation off-diagonals | COR |
| THETASE | Theta standard errors | VAR(R/2,A) |
| OMEGASE_D | Omega diagonal standard errors | VAR(R/2,A) |
| OMEGASE_O | Omega off-diagonal standard errors | COV(R/2,A) |
| SIGMASE_D | Sigma diagonal standard errors | VAR(R/2,A) |
| SIGMASE_O | Sigma off-diagonals standard errors | COV(R/2,A) |
| OMEGACSE_D | Omega correlation diagonal standard errors | VAR(R/2,A) |
| OMEGACSE_O | Omega correlation off-diagonal standard errors | COV(R/2,A) |
| SIGMACSE_D | Sigma correlation diagonal standard errors | VAR(R/2,A) |
| SIGMACSE_O | Sigma correlation off-diagonal standard errors | COV(R/2,A) |
| THETANP | Nonparametric THETAs | GENERAL |
| EXNPETA | EX non-parametric etas | GENERAL |
| COVNPETA_D | Covariance of nonparametric etas, diagonals | DIAG |
| COVNPETA_O | Covariance of nonparametric etas, off-diagonals | OFFDIAG |
| OMEGANP_D | Omega of nonparametric analysis diagonals | DIAG |
| OMEGANP_O | Omega of nonparametric analysis off-diagonals | OFFDIAG |
| COVNPETAC_D | Correlation of nonparametric etas, diagonals | DIAG (R/2,A) |
| COVNPETAC_O | Correlation of nonparametric etas, off-diagonals | COR |
| OMEGANPC_D | Omega correlation of nonparametric analysis diagonals | DIAG (R/2,A) |
| OMEGANPC_O | Omega correlation of nonparametric analysis off-diagonals | COR |
| COVARIANCE_D | Diagonals of variance-covariance of estimates | VAR |
| COVARIANCE_O | Off-Diagonals of variance-covariance of estimates | COV |
| CORRELATION_D | Diagonals of correlation of variance-covariance of estimates | VAR(R/2,A) |
| CORRELATION_O | Off-Diagonals of correlation of variance-covariance of estimates | COR |
| INVCOVARIANCE_D | Diagonals of inverse of variance-covariance of estimates | VAR |
| INVCOVARIANCE_O | Off-Diagonals of inverse of variance-covariance of estimates | COV |
| SMATRIX_D | Diagonals of S-MATRIX | VAR |
| SMATRIX_O | Off-diagonals of S-MATRIX | COV |
| RMATRIX_D | Diagonals of R-MATRIX | VAR |
| RMATRIX_O | Off-diagonals of R-MATRIX | COV |
Because of the versatility of selecting which items are to be compared and with what precision, the xml_compare program can be used for batch processing installation qualification procedures, in comparing NONMEM results of a test run against a reference run. All results given in the standard NONMEM output file are also reported in the XML file.
For example, you may wish to compare your results for example1 against the results given in the ..\examples directory of your NONMEM installation, run from your run directory, or a special installation qualification directory you may have set up:
Nmfe75 example1.ctl example1.res
xml_compare \nonmem7.2.0\examples\examples1.xml example1.xml example1.xtl >example1.difexample1.xtl would be a file you may have modified from xml_compare.xtl to suit your installation qualification needs. These .xtl files are listed in the ..\examples directory, and are simply replicates of xml_compare.xtl. You may change these for each example problem as needed. The file example1.dif will contain a list of differences, if any.
Available in the ..\util directory are some example batch processing installation files, that will execute example1 through example10l, then perform an installation qualification on these results files, against the ones in NONMEM’s ..\examples directory:
Call example.bat (this will take many hours)
Call iq.bat (this will take 10 minutes)The iq.bat repeatedly calls dif.bat. Remember to modify the “dir” option in iq.bat to point to the actual NONMEM installed directory. Also, modify dif.bat and iq.bat as needed for your particular environment. The iq.bat script will return a total differences count among all the example files. This is a convenient way of automating an installation qualification.
5 finedata Utility Program (NM73)
The utility program finedata in the ..\util directory will augment an NM-TRAN data file to incorporate additional, non-observation, time values spaced at regular increments so that when a table is generated, NONMEM can fill these records with predicted values, from which smooth prediction curves may be plotted.
The syntax is as follows:
finedata fineplot.ctlor you may use the re-direction syntax:
finedata <fineplot.ctlwhere ..\util\fineplot.ctl is an example control stream file with special commands for the finedata program. The fineplot.ctl example is extracted from part of example6.ctl:
$PROB RUN# example6 (from r2compl)
$INPUT C SET ID JID TIME DV=CONC DOSE=AMT RATE EVID MDV CMT
$DATA example6.csv IGNORE=C
$FINEDATA TSTART=0 TSTOP=50 NEVAL=100 AXIS=TIME(LIN) CMT=1,3
FILE=example6b.csvThe only records that finedata pays attention to is $INPUT, from which it obtains the column names, $DATA, from which it obtains the input data file, $FINEDATA, which contains instructions of how to fill in with additional fine increment time records, and $PROB by which problems are separated. All other control stream records are ignored. Thus, a way to create a control stream is to copy the first records describing the data layout from an existing NONMEM control stream file, and then adding the $FINEDATA record.
5.1 $FINEDATA Options
TSTART=start time (real number or integer) for creating incremental time records. If you specify FIRST, or do not specify a value for TSTART, then the time of the first record of the subject or occasion (see OCC below) is used, or when the time is less than that of the previous record, or when EVID=3 or EVID=4. If TSTART is not a number and is not FIRST, then it is interpreted as the column name in the original data set containing the start time. In such cases, the TSTART value of the first data record of the subject is used, or of the first data record, or upon occasion change (if OCC= was given), or if EVID=3 or 4, or after a re-initialization of time (indicated by the time in the data record being less than that of the previous record). Thus, TSTART could differ according to instance. The same holds true for TSTOP, TDELTA, or NEVAL (see below) if they are obtained from the data file.
OCC=name of occasion column. This is optional, and will restart the time incrementing when the occasion changes, in addition to the other conditions listed above.
NEVAL=number of incremental time records per subject (integer, or truncated if real). If not a number, then column name in the data set containing NEVAL value. If NEVAL=-1, then you wish to interpolate covariate values in the original data set, but not add any additional records.
TDELTA: Alternative to entering NEVAL, the increment in time may be entered. If not a number, then the column name in the original data set containing the TDELTA is used.
TSTOP=stop time (real number or integer) for creating incremental time records. IF TSTOP is not specified, then default is LAST, and the last record of the subject or occasion or time section is used. If TSTOP is not a number and is not LAST, then it is assumed be the column name in the original data set containing the stop time.
FILE=output data file name, to contain original data records interspersed with incremental time records.
AXIS=Name of column containing times, usually TIME. Optionally, designate (LIN) or (LOG) in parenthesis, to indicate linear or geometric time incrementing.
- If
LIN: additive time increment=(tstop-tstart)/(neval+1) - If
LOG: multiplicative time increment=(tstop/tstart)**(1/(neval+1))
DELIM=delimiter of output data file, if it is to be different from the input data file. DLEIM=S is space, DELIM=t is tab.
ITEM=number list of values for data item ITEM for which there is to be a record at each time increment. This can be done for a series of data items. For example, if you enter
$FINEDATA CMT=1,3 EVID=2,2then two records per time point are inserted, one with CMT=1, EVID=2, and the other with CMT=3, EVID=2.
Or,
$FINEDATA CMT=1,1,3,3 EVID=0,2,0,2Inserts four records per time point, with the following CMT, EVID values, in the order specified:
CMT EVID
10
12
30
32MISSING=comma-delimited-list of missing symbols. By default a period (.) and space (s) are considered missing values. Values such as 0 or -99 may be present in the data as symbols for missing values. They may be described with MISSING=0 or MISSING=-99. During interpolation, missing values will be skipped, and only records with non-missing values will be used for interpolation.
If NEVAL/=-1, only the inserted records will have filled in interpolated values, and the original records will remain untouched. When NEVAL=-1, then original records will be filled in for the specified items, but no inserted records will be added. Thus, filling missing values in original records is done as a separate action from inserting records. They may not be done simultaneously in finedata with a single $PROB, but these two actions can be accomplished by two sequential $PROB records. See finetest7.ctl to first fill in original records with interpolated values, followed by using the resulting data file as the input for the next $PROB, in which additional records are inserted:
$PROB RUN# example6 (from r2compl)
$INPUT C SET ID JID TIME DV=CONC DOSE=AMT RATE EVID MDV CMT WT
$DATA finetest.csv IGNORE=C
$FINEDATA NEVAL=-1 AXIS=TIME(LIN) MISSING=-99 WT=LIN
file=finetest7.csv
$PROB RUN# example6 (from r2compl)
$INPUT C SET ID JID TIME DV=CONC DOSE=AMT RATE EVID MDV CMT WT
$DATA finetest7.csv IGNORE=C
$FINEDATA tstart=0 TSTOP=50 NEVAL=250 AXIS=TIME(LIN) CMT=1,3 WT=LIN,PREV MISSING=-99
file=finetest7a.csv5.2 Interpolation schemes
A scheme to determine how to supply values to various data items for these inserted records may also be given. For example, to specify that the value of the next original record should be used to supply the value for WT in the inserted record:
$FINEDATA WT=NEXTThe following values may be given:
NEXT: When inserting records between two consecutive original records of time t1 (PREV) and t2 (NEXT), the PREDPP’s default of using the covariate value of the t2 (NEXT) record is used for the inserted records. NEXT is the default.
PREV: When inserting records between two consecutive original records of time t1 (PREV) and t2 (NEXT), the covariate value of the t1 (PREV) record is used for the inserted records. (LAST may be coded instead of PREV, to be consistent with the options of the $BIND record. Note that the $BIND record is not used by finedata.)
LIN, or LINLIN: A covariate-linear, time-linear interpolation is used for the covariate value for the inserted records. LINT or LINLINT (T for truncate) produces truncated integer values, LINR or LINLINR (R for round) produces values rounded to the nearest integer.
LOG, or LOGLIN: A covariate-logarithmic, time-linear interpolation is used for the covariate value for the inserted records. A T or R suffix results in truncated or rounded integer values, respectively.
LINLOG: A covariate-linear, time-logarithmic interpolation is used for the covariate value for the inserted records. A T or R suffix results in truncated or rounded integer values, respectively.
LOGLOG: A covariate-logarithmic, time-logarithmic interpolation is used for the covariate value for the inserted records. A T or R suffix results in truncated or rounded integer values, respectively.
Another example:
$FINEDATA CMT=3,3 EVID=NEXT,2indicating to create two inserted records for a given fine time point. For the first inserted record, CMT=3, and EVID of the next original record. For the second inserted record, CMT=3 and EVID=2.
Inserted records will be given the following values by default (unless over-ridden by a data item specification, such as $FINEDATA EVID=2):
DV=.
EVID=0
MDV=1Times may be entered as numerical values, or in hh:mm:ss format. Data sets with DATE/TIME records may also be processed (but then TSTART and TSTOP must be in numerical hours or hh:mm:ss format).
5.3 Using augmented data
Once finedata produces the augmented data file, in this example example6b.csv, then, a suitable NM-TRAN control stream file that would take advantage of these augmented records would be (taken from example6b.ctl in the ..\util directory):
$PROB RUN# example6 (from r2compl)
$INPUT C SET ID JID TIME DV=CONC DOSE=AMT RATE EVID MDV CMT
$DATA example6b.csv IGNORE=C
$SUBROUTINES ADVAN13 TRANS1 TOL=4
$MODEL NCOMPARTMENTS=3
$PK
...
$DES
...
$ERROR
CALLFL=0
ETYPE=1
IF(CMT.NE.1) ETYPE=0
IPRED=F
Y = F + F*ETYPE*EPS(1) + F*(1.0-ETYPE)*EPS(2)
…
$EST METHOD=ITS INTERACTION SIGL=4 NITER=25 PRINT=1 FILE=example6.ext NOABORT
$TABLE ID TIME CONC IPRED CMT MDV EVID NOAPPEND NOPRINT FILE=example6b.fin
FORMAT=,1PE12.5 ONEHEADEROf importance here is the $TABLE record. The file example6b.fin is generated by NONMEM, providing individual predicted values for each incremental time because of their presence in the input data file example6b.csv. Because incremental time records have MDV=1, there will be no impact on the estimation results. The table structure and contents of example6b.fin is suitable for importing into plotting programs, which can present smooth prediction curves (choose connect-line and no symbol) superimposed on observed data (choose with symbol, and no connect-line).
Although the added MDV=1 fine-date lines do not impact the estimation results (except where NONMEM may utilize time-changing covariates, and pick up a covariate value from these new records), they can increase estimation time. It may therefore be of advantage to perform the estimation using the original data file, followed by table generation using the enhanced data file. The FNLETA=2 setting comes in handy for this purpose:
$PROB RUN# example6 (from r2compl)
$INPUT C SET ID JID TIME DV=CONC DOSE=AMT RATE EVID MDV CMT
$DATA example6.csv IGNORE=C ; original data file used
$SUBROUTINES ADVAN13 TRANS1 TOL=4
$MODEL NCOMPARTMENTS=3
$PK
...
$DES
...
$ERROR
CALLFL=0
ETYPE=1
IF(CMT.NE.1) ETYPE=0
IPRED=F
Y = F + F*ETYPE*EPS(1) + F*(1.0-ETYPE)*EPS(2)
…
$EST METHOD=ITS INTERACTION SIGL=4 NITER=25 PRINT=1 FILE=example6.ext NOABORT
MSFO=example6.msf ATOL=4 FNLETA=0
$PROB RUN# example6 (from r2compl)
$INPUT C SET ID JID TIME DV=CONC DOSE=AMT RATE EVID MDV CMT
$DATA example6b.csv IGNORE=C ; enhanced data file
$MSFI example6.msf
$EST METHOD=1 FNLETA=2 ATOL=4
; Because FNLETA=2, no estimation us actually done. The etas loaded from the MSF file
; are used without modification to compute individual model parameters.
; Since no analysis is performed, setting METHOD=1 is sufficient, regardless of
; what method was used in the earlier analysis.
; Because ATOL=4 in the previous analysis, good idea to retain this setting, to yield
; identical evaluations from the differential equation solver.
$TABLE ID TIME CONC IPRED CMT MDV EVID NOAPPEND NOPRINT FILE=example6b.fin
FORMAT=,1PE12.5 ONEHEADERAs of NM73, if an MDV is set to a value greater than or equal to 100, it is converted to that value minus 100 upon input, but will also not be used at all during estimation, only for table outputting. This option allows you to use the same enhanced data file for estimation and Table outputs, without significantly slowing down the estimation. So, the finedata control stream file would be:
$PROB RUN# example6 (from r2compl)
$INPUT C SET ID JID TIME DV=CONC DOSE=AMT RATE EVID MDV CMT
$DATA example6.csv IGNORE=C
$FINEDATA TSTART=0 TSTOP=50 NEVAL=100 AXIS=TIME(LIN)
CMT=1,3 MDV=101,101 FILE=example6b.csvIn the following example, TSTART, TSTOP, and NEVAL are obtained from columns TIMESTART, TIMESTOP, and NEVAL, respectively.
$PROB RUN# example6 (from r2compl)
$INPUT C SET ID JID TIME DV=CONC DOSE=AMT RATE EVID MDV CMT TIMESTART
TIMESTOP NEVAL
$DATA example6c.csv IGNORE=C
$FINEDATA TSTART=TIMESTART TSTOP=TIMESTOP NEVAL=NEVAL AXIS=TIME(LIN) CMT=1,3
FILE=example6d.csvMultiple data sets may be processed by one finedata control stream file, by using $PROB records to separate the problems:
$PROB
$INPUT C=DROP ID TIME CMT OBSV DV COHT EVID AMT DOSE MDV
$DATA mydata.csv IGNORE=C
$FINEDATA tstart=0 TSTOP=700 NEVAL=500 AXIS=TIME(LIN) CMT=1,4
file=mydata_fine.csv
$PROB
$INPUT C=DROP ID TIME CMT OBSV DV COHT EVID AMT DOSE MDV
$DATA mydatab.csv IGNORE=C
$FINEDATA tstart=0 TSTOP=700 NEVAL=500 AXIS=TIME(LIN) CMT=1,4
file=mydatab_fine.csvSee also fine1, infn1, infn2 in the examples section of on-line help and guide VIII on using the INFN routine and finedata utility to create interpolated values.
5.4 Dose positions and $EXTRADOSE (NM74)
As of NM74, additional records will be inserted at all end of infusion (AMT/RATE) and additional dose positions (II, ADDL), so that predicted values will appear in the tables at these positions of discontinuity. To prevent these additional records from being inserted, set FILLDOSE=0. The dose positions and end of infusion cannot take into account any dynamic model based time-lags (ALAG), fraction infused (Fn), model dependent rates (Rn), or other model dependent changes in the actual dose events. However, if you are using RATE=-2, so that you prefer to specify duration, and the duration value is static, then you may incorporate a data item that contains the static duration, and inform finedata of it, so it can calculate where end of infusion will be for a RATE=-2 condition. Furthermore, if your fraction infused/absorbed (Fn) is static, you can add a data item that contains the static fraction infused, so that finedata can calculate end of infusion. Keep in mind that the infusion period calculated by NONMEM is (AMT*Fn)/RATE, so finedata will use the same calculation. To specify these data items for finedata to use:
LAG=lag time data itemDUR=duration data itemBIO=fraction infused data item
Here is an example for giving additional information for dosing positions (multidose.ftl):
$PROB RUN# example6 (from r2compl)
$INPUT C SET ID OC2 TIME DV=CONC DOSE=AMT RATE EVID MDV CMT ADDL II PER FF
LAGTIME
$DATA multidose.csv IGNORE=C
$FINEDATA tstart=0 TSTOP=500 NEVAL=50 AXIS=TIME(LIN) CMT=1
file=multidoseb.csv MDV=101 DUR=PER BIO=FF LAG=LAGTIMEA new option available as of NM74 is the EXTRADOSE record for finedata. EXTRADOSE adds additional non-observation dose records to those already existing, but it allows you to vary compartment number, or EVID, etc. For example, in example delayed.ftl:
$PROB Time delay problem
$INPUT ID TIME DV AMT RATE CMT EVID MDV
$DATA delayed_pre.csv IGNORE=C
$EXTRADOSE TSTART=0.0 TSTOP=50.0 AXIS=TIME CMT=3
file=delayed.csvall existing non-observation records (those for which EVID<>0) will be supplemented with additional records in which all items are identical, except the CMT will be 3, for the additional record. All other items in the record will remain unchanged from the original record for the four new records. Or,
$PROB Time delay problem
$INPUT ID TIME DV AMT RATE CMT EVID MDV
$DATA delayed.csv IGNORE=C
$EXTRADOSE TSTART=0.0 TSTOP=50.0 AXIS=TIME CMT=3,4,5,6 AMT=1,1,1,1
RATE=0,0,0,0 file=delayedb.csvwill set AMT to 1 and RATE to 0 for four new dose records with CMT numbers 3,4,5,6, consecutively. The EXTRADOSE has particular value when adding time-delay compartments to a model, to deal with time-delay problems via method of steps.
6 doexpand Utility Program for Expanding Repetitive Code (NM74)
The utility program doexpand in the ..\util directory will expand an NM-TRAN control stream file that has been annotated with DOE (which stands for DO expand) and ENDDOE (which stands for ENDDO expand) directives. This is useful for repetitive code statements that differ only by increments of indices. For example, consider the following example for replicating derivative equations that differ only by their indices to DADT() and A(), when dealing with time delay problems (although see ddexpand utility next section for a complete delay equation system process):
$DES
DOE (I=1,3,2) (J=2,4,2) (K=11,14,3)
DADT(I) = -KEL[K]*A(I)
DADT(J) = K0*(1.0 + SMAX*A(I)/(SC50*V + A(I))) - K1*A(J) &
+ KS*A(J)*SIN(2.0d+00*PI*T/TDELAY)*TDELAY/(2.0d+00*PI)
ENDDOEThis is expanded by the doexpand utility to:
DADT(1) = -KEL11*A(1)
DADT(2) = K0*(1.0 + SMAX*A(1)/(SC50*V + A(1))) - K1*A(2) &
+ KS*A(2)*SIN(2.0d+00*PI*T/TDELAY)*TDELAY/(2.0d+00*PI)
DADT(3) = -KEL14*A(3)
DADT(4) = K0*(1.0 + SMAX*A(1)/(SC50*V + A(3))) - K1*A(4) &
+ KS*A(4)*SIN(2.0d+00*PI*T/TDELAY)*TDELAY/(2.0d+00*PI)The DOE statement has index statements, of the form (VAR=START,END,STEP). In the above example, the variable I is replaced with starting value 1, ending value 3, by steps of 2. Index variable J is replaced with starting value 2, ending value 4, steps of 2, etc. The number of items in each list for each variable need not be the same. Thus, records between DOE and ENDOE are done twice, first for I=1, J=2, then for I=3, J=4. You may enter up to 200 index lists for a given DOE record. Furthermore, surround substitution variable with brackets [] when the variable is part of a variable name rather than an index, as shown in the above example for [K]. You may also nest DOE/ENDOE blocks. Negative steps are allowed, but END must then be less than START. You may use & to continue on the next line if needed. If END is not specified, END is assumed to be equal to START. IF STEP is not specified, it is assumed to be 1.
The command syntax is as follows:
doexpand original.ctl new.ctlor
doexpand <original.ctl >new.ctlIt is new.ctl that is then to be submitted to NMTRAN, via the nmfe7* script. The do expand method may also be used when providing a series of transit compartment differential equations that differ only by their DADT() and A() indexing.
An alternative specification of indexing is to use key word REPS for number of repetitions, and STEP for step size, and then only the starting index need be specified for each index variable. For example,
DOE (I=1,3,2) (J=2,4,2)could be written as
DOE (REPS=2) (STEP=2) (I=1) (J=2)See the delayed* files in the example directory for a simple working example. The batch file delayed.bat contains the following commands:
doexpand <delayed_pre.ctl >delayed.ctl
finedata <delayed.ftl
call nmfe75 delayed.ctl delayed.res –prdefaultIt first executes finedata on the instruction file delayed.ftl, to add an extra dose for the time delay compartment number 3 in data file delayed_pre.csv, produce the new data file delayed.csv. Next, doexpand acts on delayed_pre.ctl containing DOE/ENDDO records, to produce delayed.ctl. Then, the problem is executed using nonmem.
Yet another syntax that could be used is [I-1], [I+2], etc., to incorporate simple offsets of an index (up to +/-100), rather than having to specify additional indices for each offset. For example, in the previouse example, J is equivalent to I+1, so J can be replaced with [I+1]:
$DES
DOE (I=1,3,2) (K=11,14,3)
DADT(I) = -KEL[K]*A(I)
DADT([I+1]) = K0*(1.0 + SMAX*A(I)/(SC50*V + A(I))) - K1*A([I+1]) &
+ KS*A([I+1])*SIN(2.0d+00*PI*T/TDELAY)*TDELAY/(2.0d+00*PI)
ENDDOEAnother use for DOEXPAND is to prepare a partial differential equations (PDE) problem into a pure ordinary differential equations (ODE) problem. For example, diffusion type components of the form
\[\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2}\]
with a second independent integration variable x is performed numerically by finite difference expanding the right hand term. Suppose you want to expand it so that there are 20 steps (the more steps, the greater the accuracy), and the starting compartment at the left boundary is 5, and the ending compartment is at the right side of the “diffusion membrane” is modeled as compartment 25. Compartments 6 through 24 will be the intermediary diffusion compartments
…
N=20
kdel=N/L ; L is "width" of membrane (or some other appropriate meaning)
kdel2=kdel*kdel
DADT(5)={user model code}+kdel2*(A(6)-A(5))
DOE (I=6,24)
DADT(I) = kdel2*(A([I+1])-2.0*A(I)+A([I-1]))
ENDDOE
DADT(25)=kdel2*(A(24)-A(25))+{user model code}
…Notice that the equation template is simply the finite difference approximation for second derivative assessment with respect to x. This accounts for bi-directional movement of material across the diffusion barrier, where
A(I-1)-A(I)Accounts for forward movement, and
A(I+1)-A(i)Accounts for backward movement. Together they are
A(I+1)-2*A(I)+A(I-1)consistent with the code portion above.
Notice the similarity of diffusion with transit compartments, but with transit compartments, there is only forward movement:
A(I-1)-A(I)Which would be expanded as:
…
N=20
kdel=N/L
DADT(5)={user model code}-kdel*A(5)
DOE (I=6,24)
DADT(I) = kdel*(A([I-1])-A(I))
ENDDOE
DADT(25)=kdel*A(24)+{user model code}
…If the transit compartment were conceived as a one-directional transfer across some physiological region of some width L, then as number of divisions N increases, the transit compartment paradigm approaches
\[\frac{\partial u}{\partial t} = -\frac{\partial u}{\partial x}\]
since the A(I-1)-A(I) term represents a backward finite difference approximation of (negative) first derivative assessment.
Because there can be a variety of PDE expansion algorithms, one may wish to retain a repository of files with these expansions. For example, file diffusion may contain the line
DADT([I])=KDEL*( A([I-1]) - 2.0*A([I]) + A([I+1]) )Which may be “include”d into a script:
DOE (I=6,24)
include diffusion
ENDDOEThe include file may contain more than one line (it may even reference include files, up to 20 nested deep). Only include statements within a DOE..ENDDOE loop will be opened and expanded by DOEXPAND.
6.1 The DOPDE Method of Modeling PDE’s
One can insert code that DOEXPAND can then expand to create more elaborate and systematic PDE systems, to incorporate one or more dimensions of flow and diffusion spaces. The ideas are taken from Schiesser (Schiesser 1991).
Here are some commands:
DOPDE This starts the pattern setup
ENDPDE This ends the pattern setup
SIZE:n
Sets the number of compartments to describe the PDE process. Typically, this is L/dx, where L is the total length, and dx is the width that each compartment represents.
START:m Then m+1 is the starting compartment for the PDE expansion.
MVAR:name name is the variable name of the rate constant of transfer. So, if it is KDEL, then KDEL should be defined in the control stream before the DOPDE..ENDPDE segment. Default is KDEL.
MVARX:name name is the variable name of the rate constant of transfer for a derivative reference. So, if it is KDELX, then KDELX should be defined in the control stream before the DOPDE..ENDPDE segment. Default is KDELX.
OFFSET:n n is the offset compartment to the pattern. (This is explained as other commands are described).
m:a,b,c,d,…
A list of index values for the mth PDE compartment’s (that is, compartment START+m) coefficients.
R:a,b,c,d,… A list of index values for the repeat compartments’ coefficients.
NMj:a,b,c,d,e… The list of indices for the N minus jth compartment
N:a,b,c,d,e… The list of indices for the Nth compartment.
UX:a,b The a is the coefficient for du/dt coefficient (to which MVARX associates as a rate constant) of the first compartment. Set to 0 if not needed. The b is the coefficient for du/dt coefficient (to which MVARX associates as a rate constant) of the nth compartment. Set to 0 if not needed.
PATTERN:dss044_00
A pattern file, in this case, dss044_00, is “included”, which contains OFFSET, 1, 2, NM2, NM1, N, R, UX descriptions for a particular purpose. In the case of dss044_00, this contains the following pattern:
OFFSET:-5
1: 0,0,0,0,0,45,-154,214,-156,61,-10
2: 0,0,0,0,10.0,-15,-4,14,-6,1
R: 0,0,0,-1,16,-30,16,-1
NM1:0,1,-6,14,-4,-15,10
N:-10,61,-156,214,-154,45This is the pattern for coefficients to A() states for the first DADT() term, second DADT() term, repeat DADT() terms (3rd through N-2), the N-1th (NM1) DADT() term, and Nth DADT() term, in the PDE expansion. This pattern matches what is done in the dss044.for file by Schiesser:
SUBROUTINE DSS044(XL,XU,N,U,UX,UXX,NL,NU)
DX=(XU-XL)/DFLOAT(N-1)
R12DXS=1./(12.0D0*DX**2)
C... WITHOUT UX (EQUATION (53))
IF(NL.EQ.1)THEN
UXX(1)=R12DXS*
1 ( 45.0D0*U(1)
2 -154.0D0*U(2)
3 +214.0D0*U(3)
4 -156.0D0*U(4)
5 +61.0D0*U(5)
6 -10.0D0*U(6))
C...
C... WITH UX (EQUATION (36))
ELSE IF(NL.EQ.2)THEN
UXX(1)=R12DXS*
1 (-415.0D0/6.0D0*U(1)
2 +96.0D0*U(2)
3 -36.0D0*U(3)
4 +32.0D0/3.0D0*U(4)
5 -3.0D0/2.0D0*U(5)
6 -50.0D0*UX(1)*DX)
END IF
C...
C... UXX AT THE RIGHT BOUNDARY
C...
C... WITHOUT UX (EQUATION (54))
IF(NU.EQ.1)THEN
UXX(N)=R12DXS*
1 ( 45.0D0*U(N )
2 -154.0D0*U(N-1)
3 +214.0D0*U(N-2)
4 -156.0D0*U(N-3)
5 +61.0D0*U(N-4)
6 -10.0D0*U(N-5))
C...
C... WITH UX (EQUATION (37))
ELSE IF(NU.EQ.2)THEN
UXX(N)=R12DXS*
1 (-415.0D0/6.0D0*U(N )
2 +96.0D0*U(N-1)
3 -36.0D0*U(N-2)
4 +32.0D0/3.0D0*U(N-3)
5 -3.0D0/2.0D0*U(N-4)
6 +50.0D0*UX(N )*DX)
END IF
UXX(2)=R12DXS*
1 ( 10.0D0*U(1)
2 -15.0D0*U(2)
3 -4.0D0*U(3)
4 +14.0D0*U(4)
5 -6.0D0*U(5)
6 +1.0D0*U(6))
UXX(N-1)=R12DXS*
1 ( 10.0D0*U(N )
2 -15.0D0*U(N-1)
3 -4.0D0*U(N-2)
4 +14.0D0*U(N-3)
5 -6.0D0*U(N-4)
6 +1.0D0*U(N-5))
DO 1 I=3,N-2
UXX(I)=R12DXS*
1 ( -1.0D0*U(I-2)
2 +16.0D0*U(I-1)
3 -30.0D0*U(I )
4 +16.0D0*U(I+1)
5 -1.0D0*U(I+2))
1 CONTINUE
RETURN
ENDSo let’s consider the pattern for the first DADT() equation in the PDE pattern:
1: 0,0,0,0,0,45,-154,214,-156,61,-10Let us suppose that START=0, then the first DADT() is one after that, DADT(1). Because OFFSET=-5, the first coefficient is for 5 compartments before 1, A(1-5+0). Since this compartment should not exist for the pattern, its coefficient is 0. Not until we get to the 6th coefficient, belonging to compartment A(1-5+5)=A(1), should the pattern for DADT(1) begin, and it is the value 45, just as in the dss044.for source code. The coefficient to A(2) for the DADT(1) equation is then -154, etc. Thus, DADT(1) with its linear combination of A() states, will be constructed according to dss044.for (using the without UX(1) option). The 2: line contains coefficients appropriate for DADT(2).
2: 0,0,0,0,10.0,-15,-4,14,-6,1The first coefficient in the list belongs to A(2-5+0), which of course should be 0. Not until the 5th coefficient, A(2-5+4)=A(1), with value of 10, should the A() state pattern for DADT(2) begin, following the coefficients for UXX() in dss044.for.
The following are coefficients for N-1 and Nth derivatives respectively (where nth DADT() would be that of SIZE+START):
NM1:0,1,-6,14,-4,-15,10
N:-10,61,-156,214,-154,45The reason for OFFSET needing to be -5 is to accommodate the six non-zero coefficients required for describing the nth DADT (for A() states N-5, N-4, N-3, N-2, N-1, and N).
The coefficients to the internal derivatives (so, derivatives 3 through N-2) are described with the R: record
R: 0,0,0,-1,16,-30,16,-1These are the coefficients in the DO LOOP in dss044.for.
If you want to use the pattern in which UX(1) and UX(N) are used (condition NL=1, and NU=1 in the dss044.for code), then this would be as shown in dss044_11:
OFFSET:-5
UX:-50,50.0
1: 0,0,0,0,0,69.1666666666667,96.0,-36.0,10.66666666667,-1.5
2: 0,0,0,0,10.0,-15,-4,14,-6,1
R: 0,0,0,-1,16,-30,16,-1
NM1:0,1,-6,14,-4,-15,10
N:0,-1.5,10.66666666667,-36.0,96.0,69.1666666666667Notice the coefficients for the UX(1) and UX(N) term are listed as first and second coefficients listed in the UX: record, respectively.
The various PDE patterns from Schiesser (Schiesser 1991) are in the files in ..\util\pde directory, given as dss020, dss036, etc. Please read the comment in the respect .for files and reference (Schiesser 1991) to understand the basis of each of these patterns.
6.2 Combination Patterns
Let’s consider some combination patterns. Consider the file flowxdiffy.txt, so named because it describes flow (transit compartment transfer) in the x direction, and no cross transiting of diffusion along the y (or radial direction if a tube), but just parallel flow:
; flow through a tube (or blood vessel), in x direction, with no difussion
; in y (or radial direction from center). In x direction, there is
; a pre-compartment 1, and postcompartment 52.
DOPDE
; X direction, non-boundaries, so connects with compartment 1, ends
; with compartment 2+10*5=52
; Classic diffusion dudt=d^2u/dx^2
SIZE:10
START:1
MVAR:KR
PATTERN:dss012nb ; transit compartment pattern
;y DIRECTION
; Classic diffusion dudt=d^2u/dx^2, boundaried
; No pattern, but there is number of compartments (SIZE) along
; y direction to be defined
SIZE:5
ENDPDETo run this problem.
Doexpand flowxnodiffy.txt flowxnodiffy.outAnd you will see the resulting code in flowxnodiffy.out, which could be cut and pasted into the $DES section of a control stream.
Another example is flowxdiffyz.txt, which has transit compartment type flow in the x direction, and diffusions along the y and z axes:
; flow through a tube, in x direction, with diffusion
; in y and z directions also. In x direction, there is
; a pre-compartment 1, and postcompartment 252.
; In y and z direction, there are boundaries.
DOPDE
; X direction, non-boundaries, so connects with compartment 1, ends
; with compartment 2+10*5*5=252
; Classic diffusion dudt=d^2u/dx^2
SIZE:10
START:1
MVAR:KR
PATTERN:dss012nb ; transit compartment pattern
;y DIRECTION
; Classic diffusion dudt=d^2u/dx^2, boundaried
SIZE:5
PATTERN:dss042_00
;z DIRECTION
; Classic diffusion dudt=d^2u/dx^2, boundaried
SIZE:5
PATTERN:dss042_00
ENDPDEFinally, an example of diffusion in the x, y, and z directions is diffxyz.txt:
; diffusion through a tube, in x direction, with diffusion
; in y and z directions also. In x direction, there is
; a pre-compartment 1, and postcompartment 252.
; In y and z direction, there are boundaries.
DOPDE
; X direction, non-boundaries, so connects with compartment 1, ends
; with compartment 2+10*5*5=252
; Classic diffusion dudt=d^2u/dx^2
SIZE:10
START:1
MVAR:KR
PATTERN:dss042nb
;y DIRECTION
; Classic diffusion dudt=d^2u/dx^2, boundaried
SIZE:5
PATTERN:dss042_00
;z DIRECTION
; Classic diffusion dudt=d^2u/dx^2, boundaried
SIZE:5
PATTERN:dss042_00
ENDPDEWith the dss* pattern files, one can construct an elaborate combination of diffusion/flow patterns with any number of dimensions which can be inserted into the NONMEM control stream. As you can see, however, PDE equations can add many compartments to a problem, and therefore solving PDE’s can potentially be very time consuming.
7 ddexpand Utility Program for Modeling Discrete Time Delays by method of Steps (NM74)
The utility program ddexpand in the ..\util directory will expand an NM-TRAN-template control stream to propagate and incorporate delay differential equations (dde) by method of steps (MOS). The template should contain the following type of information, as exampled here with ..\examples\dde\ra2.dde. The comments explain the different components of the dde system.
; Pre-Control stream template ra2.dde used by ddexpand program to
; form functional NMTRAN control stream ra2.ctl
$PROBLEM DDE Problem: RA2
; the data file should have only DOSE input records pertaining to base equations.
; Also, the CMT must be a data item
; The ddexpand program, using finedata's EXTRADOSE facility, will add doses for
; additional compartments and call the new data file ra2_dde.csv
$INPUT ID TIME AMT RATE CMT EVID MDV DV
$DATA ra2.csv IGNORE=C
$SUBROUTINES ADVAN13 TRANS1 TOL=12
$MODEL NCOMPARTMENTS=8 ; Number of compartments will be adjusted by ddexpand
$PK
CEVID = EVID
IF (CMT /= 1) CEVID = 1
K10 = THETA(1) + ETA(1)
K12 = THETA(2) + ETA(2)
K21 = THETA(3) + ETA(3)
V1 = THETA(4) + ETA(4)
K1 = THETA(5) + ETA(5)
K2 = THETA(6) + ETA(6)
K4 = THETA(7) + ETA(7)
K5 = THETA(8) + ETA(8)
SIG1 = THETA(9) + ETA(9)
SIG2 = THETA(10) + ETA(10)
SIG3 = THETA(11) + ETA(11)
; TAU1, TAU2, TAU3, etc. are time delays. This sample has one time delay, TAU1
TAU1 = THETA(12) + ETA(12)
I0 = THETA(13) + ETA(13)
K3 = 5.0
AA = 1.0
BB = 0.5
; Set initial conditions for Base equations
A_0(1) = AA
A_0(2) = I0
A_0(6) = AA
A_0(7) = I0
; Any propagations of initial conditions and ALAG's will be placed here by ddexpand program.
$DES
; AD_x_y is the State value of A(x) delayed for time TAUy.
; These are used in the differential equations later on.
; AP_x_y is the State value of A(x) in the past, for time delay TAUy.
; That is, when T<Tauy, the AP_x_y defines A(x)
; For every AD_x_y used in the differential equations, there must be an AP_X_Y defined.
; If past is constant, then it can be as simple as AP_x_y=Initial condition constant
; (same value as A_0(x) is set to).
; Make sure AP_x_y is a function of T: do not use T-TAUy,
; as this will be done by the ddexpand program.
AP_1_1 = AA * EXP(BB * T)
AP_6_1 = AA * EXP(BB * T)
; BASE EQUATIONS ENTERED BY USER.
; Note use of AD_1_1 and AD_6_1, which warrants an expansion.
DADT(1) = K3 - (K1 / K2) * (1.0 - EXP(-K2 * T)) * A(1)
DADT(2) = K4 * A(1) - K4 * AD_1_1
DADT(3) = K4 * AD_1_1 - K5 * A(3)
CC = A(4) / V1
EFFECT = CC * (SIG1 * EXP(-SIG2 * CC) + SIG3)
DADT(4) = -K10 * A(4) - K12 * A(4) + K21 * A(5)
DADT(5) = K12 * A(4) - K21 * A(5)
DADT(6) = K3 - EFFECT * A(6) - K1 / K2 * (1.0 - EXP(-K2 * T)) * A(6)
DADT(7) = K4 * A(6) - K4 * AD_6_1
DADT(8) = K4 * AD_6_1 - K5 * A(8)
; Any delay equations necessary are placed here by the ddexpand program.
$ERROR
A1 = A(1)
A2 = A(2)
A3 = A(3)
A4 = A(4)
A5 = A(5)
A6 = A(6)
A7 = A(7)
A8 = A(8)
A9 = A(9)
A10 = A(10)
A11 = A(11)
A12 = A(12)
Y1 = A(2) + A(3)
Y2 = A(7) + A(8)
Y3 = A(3)
Y4 = A(8)
IF (CMT == 1) IPRED = Y1
IF (CMT == 2) IPRED = Y2
IF (CMT == 3) IPRED = Y3
IF (CMT == 4) IPRED = Y4
IF (CMT == 1) Y = IPRED * (1.0 + EPS(1))
IF (CMT == 1) Y = IPRED * (1.0 + EPS(2))
IF (CMT == 2) Y = IPRED * (1.0 + EPS(3))
IF (CMT == 4) Y = IPRED * (1.0 + EPS(4))
$THETA
0.32544 ; 1: K10
2.6496 ; 2: K12
2.5944 ; 3: K21
0.02645 ; 4: V
0.456 ; 5: K1
0.169 ; 6: K2
0.185 ; 7: K4
0.031 ; 8: K5
0.328 ; 9: SIG1
0.328 ; 10: SIG2
0.025 ; 11: SIG3
10.6 ; 12: TAU1
2.83 ; 13: I0
$OMEGA (0.0 FIXED)X13
$SIGMA (0.04)X4
$SIMULATION (567811 NORMAL) (2933012 UNIFORM) ONLYSIMULATION SUBPROBLEMS=1
$TABLE TIME Y1 Y2 Y3 Y4 EXCLUDE_BY CEVID NOAPPEND NOPRINT FILE=ra2.tabThus, the user need only populate the template file with the base equations above, then execute ddexpand as follows:
ddexpand ra2.dde ra2.ctlThe destination file ra2.ctl will be produced, and a finedata control stream file (fine.ftl) will be produced and submitted to finedata utility to add extra doses to the data file identified in the ra2.dde control stream (make sure finedata is in the execution directory or in the PATH system). The resulting nmtran ready file will be produced (ra2.ctl) for this example:
; Pre-Control stream template ra2.dde used by ddexpand program to
; form functional NMTRAN control stream ra2.ctl
$PROBLEM DDE Problem: RA2
; The data file should have only DOSE input records pertaining to base equations.
; Also, the CMT must be a data item.
; The ddexpand program, using finedata's EXTRADOSE facility, will add doses for
; additional compartments and call the new data file ra2_dde.csv
$INPUT ID TIME AMT RATE CMT EVID MDV DV
$DATA ra2_dde.csv IGNORE=C
$SUBROUTINES ADVAN13 TRANS1 TOL=12
$MODEL NCOMPARTMENTS=12 ; Number of compartments will be adjusted by ddexpand
$PK
CEVID = EVID
IF (CMT /= 1) CEVID = 1
K10 = THETA(1) + ETA(1)
K12 = THETA(2) + ETA(2)
K21 = THETA(3) + ETA(3)
V1 = THETA(4) + ETA(4)
K1 = THETA(5) + ETA(5)
K2 = THETA(6) + ETA(6)
K4 = THETA(7) + ETA(7)
K5 = THETA(8) + ETA(8)
SIG1 = THETA(9) + ETA(9)
SIG2 = THETA(10) + ETA(10)
SIG3 = THETA(11) + ETA(11)
; TAU1, TAU2, TAU3, etc. are time delays. This sample has one time delay, TAU1
TAU1 = THETA(12) + ETA(12)
I0 = THETA(13) + ETA(13)
K3 = 5.0
AA = 1.0
BB = 0.5
; Set initial conditions for Base equations
A_0(1) = AA
A_0(2) = I0
A_0(6) = AA
A_0(7) = I0
; Any propagations of initial conditions and ALAG's will be placed here by ddexpand program.
; INITIALIZING EQUATIONS FOR DDE COMPARTMENTS
A_0(9) = AA
A_0(12) = AA
TAU_1 = 1 * TAU1
MTDIFF = 1
MTIME(1) = PASTZERO + TAU_1
DTAU_1 = 0.0
IF ((TSTATE - PASTZERO) >= TAU_1) DTAU_1 = 1.0
CTAU_1 = 0.0
IF (TSTATE - PASTZERO >= TAU_1) CTAU_1 = 1.0
ALAG9 = TAU_1
ALAG10 = TAU_1
ALAG11 = TAU_1
ALAG12 = TAU_1
$DES
; AD_x_y is the State value of A(x) delayed for time TAUy.
; These are used in the differential equations later on.
; AP_x_y is the State value of A(x) in the past, for time delay TAUy.
; That is, when T<Tauy, the AP_x_y defines A(x)
; For every AD_x_y used in the differential equations,
; there must be an AP_X_Y defined.
; If past is constant, then it can be as simple as AP_x_y=Initial
; condition constant
; (same value as A_0(x) is set to).
; Make sure AP_x_y is a function of T: do not use T-TAUy,
; as this will be done by the ddexpand program.
AP_1_1 = AA * EXP(BB * (T - TAU_1))
AP_6_1 = AA * EXP(BB * (T - TAU_1))
AD_1_1 = (1.0 - DTAU_1) * AP_1_1 + DTAU_1 * A(9)
AD_6_1 = (1.0 - DTAU_1) * AP_6_1 + DTAU_1 * A(12)
; BASE EQUATIONS ENTERED BY USER. Note use of AD_1_1 and AD_6_1,
; which warrants an expansion.
DADT(1) = K3 - (K1 / K2) * (1.0 - EXP(-K2 * T)) * A(1)
DADT(2) = K4 * A(1) - K4 * AD_1_1
DADT(3) = K4 * AD_1_1 - K5 * A(3)
CC = A(4) / V1
EFFECT = CC * (SIG1 * EXP(-SIG2 * CC) + SIG3)
DADT(4) = -K10 * A(4) - K12 * A(4) + K21 * A(5)
DADT(5) = K12 * A(4) - K21 * A(5)
DADT(6) = K3 - EFFECT * A(6) - K1 / K2 * (1.0 - EXP(-K2 * T)) * A(6)
DADT(7) = K4 * A(6) - K4 * AD_6_1
DADT(8) = K4 * AD_6_1 - K5 * A(8)
;Any delay equations necessary are placed here by the ddexpand program.
DADT(9) = CTAU_1 * (K3 - (K1 / K2) * (1.0 - EXP(-K2 * (T - TAU_1))) * A(9))
CC_1 = A(10) / V1
EFFECT_1 = CC_1 * (SIG1 * EXP(-SIG2 * CC_1) + SIG3)
DADT(10) = CTAU_1 * (-K10 * A(10) - K12*A(10)+K21*A(11))
DADT(11) = CTAU_1 * (K12 * A(10) - K21 * A(11))
DADT(12) = CTAU_1 * (K3 - EFFECT_1 * A(12) - K1 / K2 * (1.0 - EXP(-K2 * (T - TAU_1))) * A(12))
; FOR FINEDATA $EXTRADOSE: CMT=1:,9,4:,10,5:,11,6:,12
$ERROR
A1 = A(1)
A2 = A(2)
A3 = A(3)
A4 = A(4)
A5 = A(5)
A6 = A(6)
A7 = A(7)
A8 = A(8)
A9 = A(9)
A10 = A(10)
A11 = A(11)
A12 = A(12)
Y1 = A(2) + A(3)
Y2 = A(7) + A(8)
Y3 = A(3)
Y4 = A(8)
IF (CMT == 1) IPRED = Y1
IF (CMT == 2) IPRED = Y2
IF (CMT == 3) IPRED = Y3
IF (CMT == 4) IPRED = Y4
IF (CMT==1) Y = IPRED * (1.0 + EPS(1))
IF (CMT==1) Y = IPRED * (1.0 + EPS(2))
IF (CMT==2) Y = IPRED * (1.0 + EPS(3))
IF (CMT==4) Y = IPRED * (1.0 + EPS(4))
$THETA
0.32544 ; 1: K10
2.6496 ; 2: K12
2.5944 ; 3: K21
0.02645 ; 4: V
0.456 ; 5: K1
0.169 ; 6: K2
0.185 ; 7: K4
0.031 ; 8: K5
0.328 ; 9: SIG1
0.328 ; 10: SIG2
0.025 ; 11: SIG3
10.6 ; 12: TAU1
2.83 ; 13: I0
$OMEGA (0.0 FIXED)X13
$SIGMA (0.04)X4
$SIMULATION (567811 NORMAL) (2933012 UNIFORM) ONLYSIMULATION SUBPROBLEMS=1
$TABLE TIME Y1 Y2 Y3 Y4 EXCLUDE_BY CEVID NOAPPEND NOPRINT FILE=ra2.tabNotice the expansion of appropriate equations. Also, notice that only the necessary equations are propagated, allowing conservation of memory and computation.
Delay differential equations and their propagation for their use with ODE’s is a studied art, and the above example is not to be considered to be a complete description of this process. The user is expected to be knowledgeable in this art in order to use the utility properly.
7.1 Recursive propagation for feedback loops
There are some models, with feedback loops for example, in which the time delay process and their equations must be propagated or recursively applied to cover a certain simulation time. In this case, the user should supply commented lines containing the minimal TAU1, TAU2, etc., that are used in the base equations, and the maximal simulation time TSTOP to be used. For the example EPO, the epo.dde file contains at the beginning the expected tine values for several TAUs, and TSTOP time:
; Pre-Control stream template
; TAU1 = 75.0
; TAU2 = 42.0
; TAU3 = 147.0
; TAU4 = 114.0
; TAU5 = 1587.0
; TAU6 = 1554.0
; TSTOP = 200.0When the ddexpand utility is applied, the original 6 differential equations will be expanded to 72 differential equations. It would be very tedious for the user to have to type in 66 additional equations that ddexpand does automatically, and with no typographical errors.
7.2 Sizing parameters
Certain sizing parameters may need to be modified to accommodate NMTRAN and NONMEM in producing and executing code (located in FSUBS.f90) of such extensive content. Some control stream statements to consider are as follows:
Second derivatives need not be produced if Laplace method is not used, or Laplace method NUMERICAL is used, so the following line may be added in the control stream:
$ABBR DERIV2=NO DERIV2=NOCOMMONFurthermore, if only simulations are produced, and there are no steady-state conditions requested, the production of first derivatives may also be prevented, as follows:
$ABBR DERIV2=NO DERIV2=NOCOMMON DERIV1=NOFor sizing, the number of compartments may need to be increased
$SIZES PC=73where PC must at least equal the number of compartments+1
The number of lag times and additional doses may need to be increased (as a rule of thumb, make it at least the number of compartments), for example:
$SIZES PAL=100The number of basic plus additional PK parameters (for ODE ADVAN’s, the number of user-defined parameters may be used) should be accounted for:
$SIZES PG=150Finally, there are a few parameters that must be increased to allow NMTRAN enough memory to build the FSUBS.f90 routine. The size requirements for these are reduced when you turn off DERIV2 and/or DERIV1 in the manner shown above with the $ABBR statements. The default values of these parameters are listed in ..\resource\SIZES.f90.
$SIZES DIMNEW=1000Related to number of intermediate variables, and can be incremented by 1000’s, until NMTRAN no longer issues an error.
$SIZES DIMTMP=500Related to the number of user-defined variables, and can be incremented by 1000’s, until NMTRAN no longer issues an error.
$SIZES LNP4=4000Related to generating the NMPRD4P common block module, and can be incremented by 1000’s starting from its default of 4000, until NMTRAN no longer issues an error.
Please see ddexpand-page-2017.pdf for further descriptions of the example problems.
8 Using the Delay differential equation Solvers with the ddexpand program for Discrete Delay Problems (nm75)
The advantage to using ddexpand to add the additional ODE’s inline in the NMTRAN control stream by method of steps is that NMTRAN can create proper partial derivative terms that will make it useful for FOCE, Laplace, and ITS. The disadvantage is that if there is a recursive process, as with the EPO model, then there can be a large number of equations added, and this may slow down an estimation process. An alternative is to use a delay differential equation (DDE) solver, such as RADAR5 (jointly developed by Ernst Hairer and Nicola Guglielmi (Guglielmi and Hairer 2001)) (ADVAN16) or DDE_SOLVER (ADVAN18) delay differential equation solvers (or RADAR5/ADVAN17 if there are equilibrium compartments). To prepare your control stream for these solvers, select ADVAN16, ADVAN17, or ADVAN18 in the $SUB record, then execute ddexpand (using example ..\examples\dde\epod.dde):
ddexpand epod.dde epod.ctlResults in the addition of the following setup equations in $DES:
; Pre-Control stream template UNCOMMENT TO USE
;DDE
;TAU1=75.0
;TAU2=42.0
;TAU3=147.0
;TAU4=114.0
;TAU5=1587.0
;TAU6=1554.0
;TSTOP=200.0
;$SIZES MAXNDRS=3 ; SUGGESTED VALUE FOR MAXNRDS.
$SIZES PC=80 DIMNEW=3000 DIMTMP=1000 PG=200
$PROB EPO
$ABBR FUNCTION DDEFUNC(VDDE,7)
$ABBR DERIV2=NO DERIV2=NOCOMMON ; DERIV1=NO
$INPUT ID AMT TIME DV EVID MDV CMT
$DATA EPO.csv IGNORE=C
$SUBROUTINES ADVAN16 TOL=9 ATOL=9
$MODEL NCOMPARTMENTS=6
$PK
KON =THETA(1)+ETA(1)
KOFF =THETA(2)
KEL =THETA(3)
KPT =THETA(4)
KTP =THETA(5)
VP =THETA(6)
KINT =THETA(7)
SMAX =THETA(8)
SC50 =THETA(9)
IMAX =THETA(10)
IC50 =THETA(11)
MCH =THETA(12)
C0 =THETA(13)
RR0 =THETA(14)
KDEG =THETA(15)
RBC0 =THETA(16)
TP1 =THETA(17)
TP2 =THETA(18)
TRET =THETA(19)
TRBC =THETA(20)
;TAUy
TAU1=TP1+TP2
TAU2=TP2
TAU3=TP1+TP2+TRET
TAU4=TP2+TRET
TAU5=TP1+TP2+TRET+TRBC
TAU6=TP2+TRET+TRBC
RET0=TRET*RBC0/(TRET+TRBC)
RBCM0=RBC0-RET0
HB0=MCH*RBC0
AT0=KPT*C0*VP/KTP
RC0=KON*RR0*C0/(KOFF+KINT)
KEPO=KEL*C0*VP+KINT*RC0*VP
KSYN=KDEG*RR0+KINT*RC0
KIN=RET0/(TRET*(1+SMAX*RC0/(SC50+RC0))**2)
; Initial conditions for Base equations
A_0(1)=C0*VP
A_0(2)=AT0
A_0(3)=RR0
A_0(4)=RC0
A_0(5)=RET0
A_0(6)=RBCM0
ETZ_1=ETA(1)
TSTOP=200.0
; INITIALIZING EQUATIONS FOR DDE COMPARTMENTS
$DES
; AD_x_y is the State value of A(x) delayed for time TAUy.
; AP_x_y is the State value of A(x) in the past, for time delay TAUy.
; PASTS
TAU_1=TAU1
AP_4_1=RC0
DTAU_1=0.0
IF((T-PASTZERO)>=TAU_1) DTAU_1=1.0
VDDE(1)=TAU_1
VDDE(2)=T
VDDE(3)=4
VDDE(4)=1
VDDE(5)=1
VDDE(6)=ETZ_1
VDDE(7)=0
AZ_4_1=DDEFUNC(VDDE)
AD_4_1=(1.0-DTAU_1)*AP_4_1+DTAU_1*AZ_4_1
TAU_14=TAU2
AP_4_2=RC0
DTAU_14=0.0
IF((T-PASTZERO)>=TAU_14) DTAU_14=1.0
VDDE(1)=TAU_14
VDDE(2)=T
VDDE(3)=4
VDDE(4)=2
VDDE(5)=1
VDDE(6)=ETZ_1
VDDE(7)=0
AZ_4_2=DDEFUNC(VDDE)
AD_4_2=(1.0-DTAU_14)*AP_4_2+DTAU_14*AZ_4_2
TAU_15=TAU3
AP_4_3=RC0
DTAU_15=0.0
IF((T-PASTZERO)>=TAU_15) DTAU_15=1.0
VDDE(1)=TAU_15
VDDE(2)=T
VDDE(3)=4
VDDE(4)=3
VDDE(5)=1
VDDE(6)=ETZ_1
VDDE(7)=0
AZ_4_3=DDEFUNC(VDDE)
AD_4_3=(1.0-DTAU_15)*AP_4_3+DTAU_15*AZ_4_3
TAU_21=TAU4
AP_4_4=RC0
DTAU_21=0.0
IF((T-PASTZERO)>=TAU_21) DTAU_21=1.0
VDDE(1)=TAU_21
VDDE(2)=T
VDDE(3)=4
VDDE(4)=4
VDDE(5)=1
VDDE(6)=ETZ_1
VDDE(7)=0
AZ_4_4=DDEFUNC(VDDE)
AD_4_4=(1.0-DTAU_21)*AP_4_4+DTAU_21*AZ_4_4
TAU_22=TAU5
AP_4_5=RC0
DTAU_22=0.0
IF((T-PASTZERO)>=TAU_22) DTAU_22=1.0
VDDE(1)=TAU_22
VDDE(2)=T
VDDE(3)=4
VDDE(4)=5
VDDE(5)=1
VDDE(6)=ETZ_1
VDDE(7)=0
AZ_4_5=DDEFUNC(VDDE)
AD_4_5=(1.0-DTAU_22)*AP_4_5+DTAU_22*AZ_4_5
TAU_23=TAU6
AP_4_6=RC0
DTAU_23=0.0
IF((T-PASTZERO)>=TAU_23) DTAU_23=1.0
VDDE(1)=TAU_23
VDDE(2)=T
VDDE(3)=4
VDDE(4)=6
VDDE(5)=1
VDDE(6)=ETZ_1
VDDE(7)=0
AZ_4_6=DDEFUNC(VDDE)
AD_4_6=(1.0-DTAU_23)*AP_4_6+DTAU_23*AZ_4_6
TAU_1=TAU1
AP_5_1=RET0
DTAU_1=0.0
IF((T-PASTZERO)>=TAU_1) DTAU_1=1.0
VDDE(1)=TAU_1
VDDE(2)=T
VDDE(3)=5
VDDE(4)=1
VDDE(5)=1
VDDE(6)=ETZ_1
VDDE(7)=0
AZ_5_1=DDEFUNC(VDDE)
AD_5_1=(1.0-DTAU_1)*AP_5_1+DTAU_1*AZ_5_1
TAU_15=TAU3
AP_5_3=RET0
DTAU_15=0.0
IF((T-PASTZERO)>=TAU_15) DTAU_15=1.0
VDDE(1)=TAU_15
VDDE(2)=T
VDDE(3)=5
VDDE(4)=3
VDDE(5)=1
VDDE(6)=ETZ_1
VDDE(7)=0
AZ_5_3=DDEFUNC(VDDE)
AD_5_3=(1.0-DTAU_15)*AP_5_3+DTAU_15*AZ_5_3
TAU_22=TAU5
AP_5_5=RET0
DTAU_22=0.0
IF((T-PASTZERO)>=TAU_22) DTAU_22=1.0
VDDE(1)=TAU_22
VDDE(2)=T
VDDE(3)=5
VDDE(4)=5
VDDE(5)=1
VDDE(6)=ETZ_1
VDDE(7)=0
AZ_5_5=DDEFUNC(VDDE)
AD_5_5=(1.0-DTAU_22)*AP_5_5+DTAU_22*AZ_5_5
TAU_1=TAU1
AP_6_1=RBCM0
DTAU_1=0.0
IF((T-PASTZERO)>=TAU_1) DTAU_1=1.0
VDDE(1)=TAU_1
VDDE(2)=T
VDDE(3)=6
VDDE(4)=1
VDDE(5)=1
VDDE(6)=ETZ_1
VDDE(7)=0
AZ_6_1=DDEFUNC(VDDE)
AD_6_1=(1.0-DTAU_1)*AP_6_1+DTAU_1*AZ_6_1
TAU_15=TAU3
AP_6_3=RBCM0
DTAU_15=0.0
IF((T-PASTZERO)>=TAU_15) DTAU_15=1.0
VDDE(1)=TAU_15
VDDE(2)=T
VDDE(3)=6
VDDE(4)=3
VDDE(5)=1
VDDE(6)=ETZ_1
VDDE(7)=0
AZ_6_3=DDEFUNC(VDDE)
AD_6_3=(1.0-DTAU_15)*AP_6_3+DTAU_15*AZ_6_3
TAU_22=TAU5
AP_6_5=RBCM0
DTAU_22=0.0
IF((T-PASTZERO)>=TAU_22) DTAU_22=1.0
VDDE(1)=TAU_22
VDDE(2)=T
VDDE(3)=6
VDDE(4)=5
VDDE(5)=1
VDDE(6)=ETZ_1
VDDE(7)=0
AZ_6_5=DDEFUNC(VDDE)
AD_6_5=(1.0-DTAU_22)*AP_6_5+DTAU_22*AZ_6_5
CC=A(1)/VP
AT=A(2)
RR=A(3)
RC=A(4)
RET=A(5)
RBCM=A(6)
X1=1+SMAX*AD_4_1/(SC50+AD_4_1)
X2=1+SMAX*AD_4_2/(SC50+AD_4_2)
X3=1+SMAX*AD_4_3/(SC50+AD_4_3)
X4=1+SMAX*AD_4_4/(SC50+AD_4_4)
X5=1+SMAX*AD_4_5/(SC50+AD_4_5)
X6=1+SMAX*AD_4_6/(SC50+AD_4_6)
X0=1+SMAX*RC0/(SC50+RC0)
I1=1-IMAX*(MCH*(AD_5_1+AD_6_1)-HB0)/(IC50+(MCH*(AD_5_1+AD_6_1)-HB0))
I3=1-IMAX*(MCH*(AD_5_3+AD_6_3)-HB0)/(IC50+(MCH*(AD_5_3+AD_6_3)-HB0))
I5=1-IMAX*(MCH*(AD_5_5+AD_6_5)-HB0)/(IC50+(MCH*(AD_5_5+AD_6_5)-HB0))
; BASE EQUATIONS.
DADT(1)=KEPO-KON*CC*VP*RR+KOFF*RC*VP-(KEL+KPT)*CC*VP+KTP*AT
DADT(2)=KPT*CC*VP-KTP*AT
DADT(3)=KSYN-KON*CC*RR+KOFF*RC-KDEG*RR
DADT(4)=KON*CC*RR-(KOFF+KINT)*RC
DADT(5)=KIN*X1*X2*I1-KIN*X3*X4*I3
DADT(6)=KIN*X3*X4*I3-KIN*X5*X6*I5And, there are no new ODE equations added to the control stream. Instead, the calls to DDEFUNC provide the time-delays for the required state variables. The DDEFUNC routine calls algorithms in the DDE solver system to create the delayed versions of state variable given by VDDE(3), for delay time given by VDDE(1). The advantage to using RADAR5 is its use in those problems that would normally create a recursive scheme if they were expanded by method of steps, potentially resulting in hundreds of additional ODE’s. With the RADAR5/ADVAN16 method, no additional ODE’s are added to the control stream.
The ..\pr\RADAR5U.f90 (for ADVAN16 and ADVAN17) and ..\pr\DDESLVU.f90 (for ADVAN18) routines offer additional control. Please read ..\guides\manrad5-v2.pdf for RADAR5, and ..\guides\ddes_f90.pdf for DDE_SOLVER for details about the parameter settings, and optional routines that may be incorporated.
For example, to use DDESLVu.f90, make a copy of it in your run directory, and rename it, such as DDESLVU2.f90. Make the modifications, and then reference it in the $SUBR record:
$SUBR … OTHER=DDESLVU2.f90Some parameters of interest regarding ADVAN16(17)/RADAR5 include those controlled by $SIZES record:
PAST_SIZE=4000 by default, which determines the number of detailed times points (and hence resolution) of the delay equation storage
MAXNRDS=PC by default, but can be set to actual number of delay compartments used, in oder to save memory.
These parameters can be set with $SIZES in the control stream file:
$SIZES MAXNRDS=2 PAST_SIZE=6000ADVAN18/DDE_SOLVER operates with allocatable arrays, and does not require pre-sizing by the user.
If you are not interested in viewing the processed control stream file with the additional records added to interface with the DDE solver system, you can directly send your control stream file to the nmfe75 script, but add the -dde option. Examples in ..\examples\dde directory, such as epodd.ctl, and rad.ctl, can be directly submitted, as follows:
nmfe75 epodd.ctl epodd.res -prdefault –ddeand the problem will run. If you wish to see the processed control stream after the execution, this is file epodd.ctl_dde in your run directory. See ..\guides\advan16_dde.pdf for more information.
8.1 PASTZERO and pre-time-zero dosing
Sometimes you may wish to have the equations transition from the past to the present other than at time T=0. In this case, set the reserved variable PASTZERO to a non-zero (including negative) value:
$PK
PASTZERO=-10.0For example, suppose you wanted to have 30 additional doses (for a total of 31) every 8 hours, followed by records sampling the decline in concentration after the last dose. Suppose also that you wanted the time of sampling begin at TIME=0, therefore the beginning of the dose would be at -240.0 hours (see example simpledii16_2):
CID TIME AMT RATE II ADDL CMT EVID MDV DV
100 -240 100 0 8 30 1 1 1 0
100 -1.00E-06 0 0 0 0 1 0 0 1
100 0 0 0 0 0 1 0 0 1
100 1 0 0 0 0 1 0 0 1
100 2 0 0 0 0 1 0 0 1
100 3 0 0 0 0 1 0 0 1
100 4 0 0 0 0 1 0 0 1
100 5 0 0 0 0 1 0 0 1
100 6 0 0 0 0 1 0 0 1When using ADVAN16, ADVAN17, and ADVAN18 the past would be any time before the first dose, so it would be -240.0 hours as well. In $PK, you specify this using the reserved variable PASTZERO:
$PK
PASTZERO=-240.08.2 Using ADVAN9, 13, 14, 15 in DDE Mode
As of nm75, the traditional ODE solvers ADVAN9, 13, 14, and 15 have a superimposed algorithm that stores state variables of past times, and interpolates them for the present time (using cubic hermite interpolation, as is done by Petzoldt and Soetaert for the R DeSolve packages), which allows you to use these in a similar manner as ADVAN16 to ADVAN18, without using method of steps. These superimposed algorithms on traditional ODE solvers are not as efficient as the specially designed solvers ADVAN16 to ADVAN18, and may require an additional increase in TOL of 1 or 2 in order to obtain the desired precision. In turn the ATOL may be increased by 1 or 2, or you may increase it even to ATOL=12-15, and this may provide significant improvement in the TOL as well, without too much increase in computation time. For example, suppose that you wish to have at least 6 digits precision in your results, so that normall you would set TOL=6 ATOL=6. For a -dde problem using ADVAN 9,13,14 or 15, you could set
$SUBROUTINES ADVAN13 TOL=8 ATOL=14to achieve a 6 digit precision. A little trial and error is warranted with one iteration of an estimation step, to ensure there is not too much of a speed sacrifice.
To ensure that ddexpand does not create a method of steps (MOS) expansion, but a dde (DDE) expansion when using ADVAN 9, 13, 14, or 15, insert ;DDE as a comment line near the top of the control stream:
;DDE
…
$SUBROUTINES ADVAN13 TRANS1 TOL=8 ATOL=8Then, execute in the same manner as you would with the dde solvers ADVAN16 to ADVAN18:
nmfe75 epodd.ctl epodd.res -prdefault –ddeSimilarly, if you wish to have ADVAN16 to 18 utilize the MOS method instead of its native DDE method, you would insert ;MOS near the top of the control stream:
;MOSOtherwise, if ;DDE or ;MOS is not placed in the control stream, the default action will be MOS for ADVAN’s 9,13,14,15, and DDE for ADVAN’s 16, 17, 18.
8.3 Sensitivity Equations with ADVAN14 / ADVAN15
When using ADVAN14 or ADVAN15 in their DDE mode, and when analytical derivatives are being used (such as during ITS, FOCE, Laplace), and there is a random eta effect (inter-subject variability) on a TAU, then it is sometimes best to uncomment the lines
NBASE=NTOTAL
NSENSE=0in CVODEU.f90 (for ADVAN14) or IDAU.f90 (for ADVAN15) to turn off the sensitivity equation facility, and add in the control stream an OTHER statement for this modified code:
$SUBROUTINES ADVAN14 TOL=6 ATOL=6 OTHER=CVODEU.F90or
$SUBROUTINES ADVAN15 TOL=6 ATOL=6 OTHER=IDAU.F90Otherwise, the derivatives with respect to eta relating to TAU may be less accurately evaluated to a sufficient degree to affect optimization.
8.4 MTDIFF / MTIME for Equilibrium Compartment ADVANs
For the equilibrium compartment ADVAN’s (ADVAN9, ADVAN15, ADVAN17), even when they are doing a non-equilibirum compartment problem, when ddexpand inserts code for the DDE mode, the start of the delay time is signaled to the solver with an addition MTIME() statement, such as:
$PK
…
MTDIFF=1
MTIME(1)=PASTZERO+TAU_1
DTAU_1=0.0
IF((TSTATE-PASTZERO)>=TAU_1) DTAU_1=1.0For the other ADVAN’s the following code is inserted in the $DES section instead:
$DES
DTAU_1=0.0
IF((T-PASTZERO)>=TAU_1) DTAU_1=1.0These decisions were made based on trial and error assessments of each of these ADVAN’s abilities to properly calculate the derivative of TAUx with respect to ETAs. If the analysis, particularly FOCE/Laplace/ITS, tends to not converge efficiently, you may wish to over-ride these default decisions with ;MTDIFF or ;NOMTDIFF comment inserted near the top of the control stream. For example, to have the MTDIFF/MTIME() code used for ADVAN18, then add
;MTDIFFAnd for quilibrium compartment ADVAN’s, to not have the MTDIFF/MTIME() code inserted, add
;NOMTDIFF8.5 Symbolic eta indices and NETAS
The ddexpand program, when it adds code for DDE solvers, needs to know the number of etas used in the model. Usually it can determine this on its own, but if you are using symbolic names for eta indices (such as ETA(CL)), in the code, then it cannot interpret this. To assist ddexpand in this, add the following commented line, preferably right after the ;DDE comment line, for example:
;DDE
;NETAS=68.6 Caveat on INITIALOFF with delay solvers
When using delay solvers (;DDE), do not use the INITIALOFF feature in $MODEL compartment descriptions, unless the compartment remains off throughout for an individual for a given occasion, or the compartment is turned on (by an input dose) at the start (usually time=0) of the numerical integration, and stays on. For example, you may use it in data sets in which an individual and occasion (EVID=1 or EVID=4) is given a non-IV dose at time 0, or the individual/occasion is never given a non-IV dose.
9 Modeling Discrete Delays on General Variables Using the ddexpand program (nm751)
Sometimes you want to impose a direct delay on a general variable, rather than a state variable A(*). Consider the following case, where it is desired that kin be delayed by an amount TAU. You can readily delay the states that KIN depends on, in this case A(1) (distrib_delay\smax_discrete.ctl), creating kindelay as the TAU delayed version of kin:
;DDE
$PROBLEM Discrete Delay on Smax Drug Responsive Input function
$ABBR DERIV2=NO DERIV2=NOCOMMON DERIV1=NO
$INPUT ID AMT TIME DV EVID MDV
$DATA smax_discrete.dat IGNORE=C
$SUBROUTINES ADVAN18 TOL=10 ATOL=10
$MODEL NCOMPARTMENTS=2
$PK
CALLFL=-2
MXSTEP=2000000000
SMAX=THETA(1)
SC50=THETA(2)
KIN0=THETA(3)
AA=THETA(4)
BB=THETA(5)
TAU1 = aa*dexp(gamln(1.0+1.0/bb))
Vmax = 2
V50 = 3
V=1
; Initial conditions
A_0(1)=0.0
A_0(2)=5.0
$DES
AP_1_1=0.0
c = A(1)/V
kin = kin0*(1.0+(Smax*c)/(SC50+c))
cdelay=AD_1_1/V
kindelay=kin0*(1.0+(Smax*cdelay)/(SC50+cdelay))
DADT(1) = -(Vmax*A(1))/(V50+A(1))
DADT(2) = kin - kindelay
$ERROR
A1=A(2)
Y1=1.0
IPRED=A(2)
Y=IPRED*(1.0+EPS(1))
$THETA
0.5
10.0
0.8
7.0
6.5
$OMEGA (0.0 fixed)x4
$SIGMA (0.0 fixed)
$SIML (122345) ONLYSIM SUBP=1
$TABLE TIME KIN IPRED NOPRINT NOAPPEND FILE=smax_discrete.tab ONEHEADER
FORMAT=S1PE20.13Alternatively, you could create an additional state variable that represents kin, such as A(3), but you would need to represent the derivative of kin with respect T, DADT(3). The following example does this (smax_discrete_kin.ctl):
;DDE
$PROBLEM Discrete Delay on Smax Drug Responsive Input function kin
$ABBR DERIV2=NO DERIV2=NOCOMMON
$INPUT ID AMT TIME DV EVID MDV
$DATA smax_discrete.dat IGNORE=C
$SUBROUTINES ADVAN16 TOL=10 ATOL=10
$MODEL NCOMPARTMENTS=3
$PK
CALLFL=-2
MXSTEP=2000000000
SMAX=THETA(1)
SC50=THETA(2)
KIN0=THETA(3)
AA=THETA(4)
BB=THETA(5)
Vmax = 2
V50 = 3
V=1
TAU1 = aa*dexp(gamln(1.0+1.0/bb))
; Initial conditions
A_0(1)=0.0
A_0(2)=5.0
; KIN at time 0 actually contains the bolus dose contribution into A(1).
; So it is not kin0, but rather;
CCC = AMT/V
A_0(3)=kin0*(1.0+(Smax*ccc)/(SC50+ccc))
$DES
c = A(1)/V
; DERIVATIVE OF A WRT T, WHICH EQUALS V TIMES DERIVATIVE OF C WRT T
DADT1=-(Vmax*A(1))/(V50+A(1))
; DERIVATIVE OF KIN WRT C
DKINC=KIN0*SMAX*SC50/(SC50+C)/(SC50+C)
; The past, is when T<0, so before the bolus dose into A(1) is given, which is:
AP_3_1=KIN0
DADT(1) = DADT1
; DERIVATIVE OF KIN WRT C TIMES DERIVATIVE OF C WRT T YIELDS DERIVATIVE OF KIN WRT T
DADT(3) = DKINC*DADT1/V
; SUBMIT STATE VARIABLE A(3), WHICH REPRESENTS KIN, FOR WEIBULL-DISTRIBUTED DELAY
DADT(2) = A(3) - AD_3_1
$ERROR
CC = A(1)/V
kin = kin0*(1.0+(Smax*cc)/(SC50+cc))
A1=A(2)
KINCALC=A(3)
Y1=1.0
IPRED=A(2)
Y=IPRED*(1.0+EPS(1))
$THETA
0.5
10.0
0.8
7.0
6.5
$OMEGA (0.0 fixed)x5
$SIGMA (0.0 fixed)
$SIML (122345) ONLYSIM SUBP=1
$TABLE TIME KIN KINCALC IPRED NOPRINT NOAPPEND FILE=smax_discrete_kin.tab
ONEHEADER FORMAT=S1PE20.13As a third option, if we use the DAE ADvans, such as ADVAN15 or ADVAN17, then we can equate the A(3) with kin, and have A(3) solved via an equilibrium status, and we will not need to manually create the derivative of kin with respect to T (smax_discrete_kin_dae.ctl):
;DDE
$PROBLEM Discrete Delay on Smax Drug Responsive Input function kin using dae
$ABBR DERIV2=NO DERIV2=NOCOMMON ; DERIV1=NO
$INPUT ID AMT TIME DV EVID MDV
$DATA smax_discrete.dat IGNORE=C
$SUBROUTINES ADVAN17 TOL=10 ATOL=10 ; OTHER=ddeslvu.f90
$MODEL NCOMPARTMENTS=3
COMP(COMP1)
COMP(COMP2)
COMP(COMP3 EQUILIBRIUM)
$PK
CALLFL=-2
MXSTEP=2000000000
SMAX=THETA(1)
SC50=THETA(2)
KIN0=THETA(3)
AA=THETA(4)
BB=THETA(5)
TAU1 = aa*dexp(gamln(1.0+1.0/bb))
Vmax = 2
V50 = 3
V=1
; Initial conditions
A_0(1)=0.0
A_0(2)=5.0
$DES
AP_3_1=kin0
DADT(1) = -(Vmax*A(1))/(V50+A(1))
DADT(2) = A(3) - AD_3_1
$AESINIT
INIT=0
$AES
c = A(1)/V
kin = kin0*(1.0+(Smax*c)/(SC50+c))
E(3)=A(3)-kin
$ERROR
A1=A(2)
Y1=1.0
IPRED=A(2)
Y=IPRED*(1.0+EPS(1))
$THETA
0.5
10.0
0.8
7.0
6.5
$OMEGA (0.0 fixed)x4
$SIGMA (0.0 fixed)
$SIML (122345) ONLYSIM SUBP=1
$TABLE TIME KIN IPRED NOPRINT NOAPPEND FILE=smax_discrete_kin_dae.tab
ONEHEADER FORMAT=S1PE20.13Although for the above simple example it was sufficiently easy to delay the source state variable A(1) to kin, there are problems in which the general variable relies on several state variables, each of which are to be delayed the same amount. Rather than evaluating for each of their delayed values individually, and then constructing the final delayed variable, in such cases, using the DAE method to equate that general variable to an equilibirium compartment allows the DDE solver to calculate that delayed variable.
These techniques will also come in handy when we use the discrete differential equation solvers to solve distributed delay problems, shown in subsequent sections.
10 Solving Distributed Delayed Output Rate Problems Using Repetition Variables RPTO/RPTI (nmVI)
There are conditions in which one would like to model a delay in a distributed manner. In one special type of such 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 as
\[ \frac{dN(t)}{dt} = k_{in}(t) - \int_0^\infty l(s) k_{in}(t - s) ds \]
The second term on the right hand side is the convolution integral, and it weights input rate delayed by a time s, kin(t-s), with the pdf function l(s), and then sums over all possible time delays, from s=0 to s=infinity. This differs from the discrete delay, where only a single delay time s=tau is considered. Mathematically, such a delay at a single time tau would be represented by the delta pdf, which equals 1 at s=tau, and 0 otherwise:
\[ l(s) = \delta(s - \tau) \]
\[ \delta(z) = 1, \quad z = 0 \]
\[ \delta(z) = 0, \quad z \neq 0 \]
So that
\[ \begin{aligned} \frac{dN(t)}{dt} &= k_{in}(t) - \int_0^\infty \delta(s - \tau) k_{in}(t - s) ds \\ &= k_{in}(t) - k_{in}(t - \tau) \end{aligned} \]
While this convolution method potentially can be used to solve discrete delay problems it would be far less efficient then using discrete delay solvers.
The state variable to a distributed delayed output rate problem type can be shown to be solved in the form of:
\[ N(t) = \int_0^\infty (1 - L(s)) k(t - s) ds \]
where L is the cumulative distribution function (CDF), or integral, of l. In guide VIII, the Repetition Parameters I example shows how a convolution for this type of problem can be performed in NONMEM, using RPTO and RPTI variables. Appendix M of the NONMEM Technical Guide further explains how the NONMEM RPTO/RPTI method can use repeated use of an ODE solver to evaluate the above equation. Here we show an example that uses the WEIBULL distribution on an input function governed by an Smax relation to the drug level in the central compartment (examples\distrib_delay\smax_weibull_rep.ctl):
$PROBLEM Weibull Distributed Delay on Smax Drug Responsive Input function
using NONMEM's repetition variables
$ABBR DERIV2=NO DERIV2=NOCOMMON ; DERIV1=NO
$INPUT ID AMT TIME DV EVID MDV
$DATA smax_weibull.dat IGNORE=C
$SUBROUTINES ADVAN13 TOL=6 ATOL=6 ; OTHER=ddeslvu.f90
$MODEL NCOMPARTMENTS=2
$INFN ;enables use of repetition feature
IF (ICALL.EQ.0) RPTO=1
$PK
CALLFL=-2
MXSTEP=2000000000
SMAX=THETA(1)
SC50=THETA(2)
KIN0=THETA(3)
AA=THETA(4)
BB=THETA(5)
Vmax = 2
V50 = 3
V=1
IF (RPTI.EQ.0) TI=TIME
IF (NEWIND.EQ.2) RPTO=-1
; Initial conditions
A_0(1)=0.0
A_0(2)=5.0
$DES
c = A(1)/V
kin = kin0*(1.0+(Smax*c)/(SC50+c))
DD=0.0
if(ti-t>=0.0) DD=EXP(-((TI-T)/BB)**AA) ; DD=1-WEIBULLCDF(TI-T)
DADT(1) = -(Vmax*A(1))/(V50+A(1))
DADT(2) = DD*(KIN-KIN0) ; A(2)=N(t)
$ERROR
A1=A(2)
Y1=1.0
IPRED=A(2)
Y=IPRED*(1.0+EPS(1))
A3=A(3)
$THETA
0.5
10.0
0.8
7.0
6.5
$OMEGA (0.0 fixed)x4
$SIGMA (0.0 fixed)
$SIML (122345) ONLYSIM SUBP=1
$TABLE TIME KIN IPRED NOPRINT NOAPPEND FILE=smax_weibull_rep.tab ONEHEADER
FORMAT=S1PE20.13The advantage to using the NONMEM repetition variables method is that the distributional convolution of kin is exactly evaluated, whereas the distributed transit or distributed discrete methods described in the next sections only approximates these convolutions. The repetition variables method does add 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. Appendix M of the NONMEM Technical Guide discusses this further, and why the disadvantage of the RPTO/RPTI method is that it cannot be used for general convolution problems, and therefore cannot be used for more general distributed delay problems, whereas distributed transit and distributed discrete delay methods, shown in the next section, can be used to solve more general convolution problems.
11 Transit Compartments for Gamma distributed Time Delays Using the ddexpand program (nm751)
There are conditions in which one would like to model a delay in a distributed manner for more general problems than the delayed output rate problem. One distribution that is popular is the gamma distribution that can be emulated by the the transit compartment method. This is done by adding several ODE’s to the NMTRAN control stream, so that a state variable A(x) is filtered through a series of NU transit compartments in which the rate constant of transfer KTR 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:
KTR_y=NUy/TAUy
This is modeled in the same way as the discrete delay, but instead of specifying AP_x_y and AD_x_y, you specify APT_x_y and ADT_x_y. The transit compartment method is equivalent to a gamma distribution when NUy is an integer. 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 (CEILING(NUy)-1) and 0.75 of transit compartment CEILING(NUy) 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 only a few percent, and is very easy to calculate. Consider the following user-created nmtran control stream (distrib_delay\logt.ctl):
;DDE
;Define as a comment the largest size for NU1 that one is likely to need.
;ddexpand will create up to 10 transfer compartments for this example.
; Specify as ;NN or ;NU
;NN1=10
$PROBLEM LOGISTIC with fractional transit compartments; turn off
second derivative assessments, sometimes even 1st derivatives if only
simulating
$ABBR DERIV2=NO DERIV2=NOCOMMON ; DERIV1=NO
$INPUT ID AMT TIME PRDV DV EVID MDV
$DATA logt.dat IGNORE=C
$SUBROUTINES ADVAN15 TOL=10 ATOL=10 ; OTHER=ddeslvu.f90
$MODEL NCOMPARTMENTS=1
$PK
CALLFL=-2
MXSTEP=2000000000
KG=THETA(1)
Y0=THETA(2)
YSS=THETA(3)
TAU1=THETA(4)
; The NU1 specified here should be no greater than the NN1 transit
compartments
; defined in the comment earlier.
NU1=2.5
; Initial conditions
A_0(1)=Y0
$DES
; AD_1_1 is the State value of A(1) delayed for time TAU1.
; AP_1_1 is the State value of A(1) in the past, for time delay TAU1.
APT_1_1=Y0
DADT(1)=KG*(1.0-ADT_1_1/YSS)*A(1)
$ERROR
A1=A(1)
Y1=1.0
IPRED=A(1)
Y=IPRED*(1.0+EPS(1))
$THETA
0.2D+00
1.0D+00
10.0D+00
5.0D+00
$OMEGA (0.0 fixed)x4
$SIGMA (0.0 fixed)
$SIML (122345) ONLYSIM SUBP=1
$TABLE TIME IPRED NOPRINT NOAPPEND FILE=logt.tab ONEHEADER FORMAT=S1PE20.13The NN1 defines the number of compartments for the expansion, and TAU1 defines the 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 a constant past), because ddexpand uses the APT_x_y entry to identify the type of delay as the simple transit compartment system. The expansion by ddexpand results in:
;DDE
;Define as a comment the largest size for NU1 that one is likely to need.
;ddexpand will create up to 10 transfer compartments for this example.
; Specify as ;NN or ;NU
;NN1=10
$PROBLEM LOGISTIC with fractional transit compartments
; turn off second derivative assessments, sometimes even 1st derivatives if
only simulating
$ABBR DERIV2=NO DERIV2=NOCOMMON ; DERIV1=NO
$INPUT ID AMT TIME PRDV DV EVID MDV
$DATA logt.dat IGNORE=C
$SUBROUTINES ADVAN15 TOL=10 ATOL=10 ; OTHER=ddeslvu.f90
$MODEL NCOMPARTMENTS=11
$PK
CALLFL=-2
MXSTEP=2000000000
KG=THETA(1)
Y0=THETA(2)
YSS=THETA(3)
TAU1=THETA(4)
; The NU1 specified here should be no greater than the NN1 transit
compartments
; defined in the comment earlier.
NU1=2.5
; Initial conditions
A_0(1)=Y0
NU_1=NU1
KTR_1=NU_1/TAU1
XNU_1=MAX(INT(NU_1)+1,0)
YNU_1=NU_1-XNU_1+1.0
$DES
; AD_1_1 is the State value of A(1) delayed for time TAU1.
; AP_1_1 is the State value of A(1) in the past, for time delay TAU1.
; DELAY SETUP FOR EQUATION SET 1
APT_1_1=Y0
TZ_GD=T-PASTZERO
TZ_G=TZ_GD+TG_DEL
T_O=0.0
IF(TZ_GD>=0.0) T_O=1.0
IF(T_O==0.0) TZ_G=1.0
IX_=1
SUM_1=0.0
IF(IX_>XNU_1) THEN
DADT(2)=0.0
ENDIF
IF(IX_==XNU_1) THEN
SUM_1=T_O*((1.0-YNU_1)*A(1)+YNU_1*(A(2)+APT_1_1))
ENDIF
IF(IX_<=XNU_1) THEN
IX_=IX_+1
DADT(2)=T_O*KTR_1*(A(1)-A(2)-APT_1_1)
ENDIF
IF(IX_>XNU_1) THEN
DADT(3)=0.0
ENDIF
IF(IX_==XNU_1) THEN
SUM_1=T_O*((1.0-YNU_1)*A(2)+YNU_1*A(3)+APT_1_1)
ENDIF
IF(IX_<=XNU_1) THEN
IX_=IX_+1
DADT(3)=T_O*KTR_1*(A(2)-A(3))
ENDIF
IF(IX_>XNU_1) THEN
DADT(4)=0.0
ENDIF
…
IF(IX_>XNU_1) THEN
DADT(11)=0.0
ENDIF
IF(IX_==XNU_1) THEN
SUM_1=T_O*((1.0-YNU_1)*A(10)+YNU_1*A(11)+APT_1_1)
ENDIF
IF(IX_<=XNU_1) THEN
IX_=IX_+1
DADT(11)=T_O*KTR_1*(A(10)-A(11))
ENDIF
ADT_1_1=(1.0-T_O)*APT_1_1+SUM_1
; DELAY EQUATIONS FOR EQUATION SET 0 (BASE EQUATIONS)
DADT(1)=KG*(1.0-ADT_1_1/YSS)*A(1)
$ERROR
A1=A(1)
Y1=1.0
IPRED=A(1)
Y=IPRED*(1.0+EPS(1))
$THETA
0.2D+00
1.0D+00
10.0D+00
5.0D+00
$OMEGA (0.0 fixed)x4
$SIGMA (0.0 fixed)
$SIML (122345) ONLYSIM SUBP=1
$TABLE TIME IPRED NOPRINT NOAPPEND FILE=logt.tab ONEHEADER FORMAT=S1PE20.1312 Creating Transit Compartments for Time Delays of any distribution, Using the ddexpand program (nm751)
A general method of modeling distributed delays (without having to use discrete delay solvers, which will be shown in the next example), in which any distribution function may be used, is the general distributed delay method, developed by Koch and Schropp (Koch and Schropp 2015). In this case 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, for mean transit time TAUy and state variable A(x), and the ddexpand utility will then make the appropriate expansions.
Consider an example for a Weibull distribution. The user-written control stream may look like this (distrib_delay\logdt.ctl):
;DDE
;NN1=28
;PRC1=0.01
;DISTRIB1=WEIBULLCDF(VG)
$PROBLEM LOGISTIC with distributed transit compartments
; turn off second derivative assessments, sometimes even 1st derivatives if
only simulating
$ABBR DERIV2=NO DERIV2=NOCOMMON ; DERIV1=NO
$INPUT ID AMT TIME PRDV DV EVID MDV
$DATA logt.dat IGNORE=C
$SUBROUTINES ADVAN15 TOL=10 ATOL=10 ; OTHER=ddeslvu.f90
$MODEL NCOMPARTMENTS=1
$PK
CALLFL=-2
MXSTEP=2000000000
KG=THETA(1)
Y0=THETA(2)
YSS=THETA(3)
TAU1=THETA(4)
NU1=2.5
; Initial conditions
A_0(1)=Y0
$DES
; Define parameters to Weibull distribution again, for $DES record
VG(2)=NU1
VG(3)=TAU1
; ADG_1_1 is the State value of A(1) delayed for time TAU1, using a general
distribution.
; APG_1_1 is the State value of A(1) in the past, for time delay TAU1.
APG_1_1=Y0
DADT(1)=KG*(1.0-ADG_1_1/YSS)*A(1)
$ERROR
A1=A(1)
Y1=1.0
IPRED=A(1)
Y=IPRED*(1.0+EPS(1))
$THETA
0.2D+00
1.0D+00
10.0D+00
5.0D+00
$OMEGA (0.0 fixed)x4
$SIGMA (0.0 fixed)
$SIML (122345) ONLYSIM SUBP=1
$TABLE TIME IPRED NOPRINT NOAPPEND FILE=logdt.tab ONEHEADER FORMAT=S1PE20.13Notice the beginning comment lines:
;NN1=27
;PRC1=0.01
;DISTRIB1=WEIBULLCDF(VG)The number of transit equations to expand is NNy, for time delay y. Next, the precision PRCy is specified. If it is 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 (Koch and Schropp 2015):
WEIBULLCDFGAMMACDFLOGNORMALCDFGOMPMAKECDF(Gompertz-Makeham)WEIBULLCCDF(combined Weibull)UNIFORMCDF
In theory, one can use any of the CDF versions of the distributions listed in section I.34 Probability Density Functions (NM742), but the CDF inverse functions (distribCDFINV) to many of these have not yet been constructed, which is required for this process.
An approximation of the convolution of the distribution with the input state variable A(x) is intended to be calculated. The integral of the CDF is to be evaluated at evenly spaced time intervals from s=0 to s= distribCDFINV(1-PRC1), spaced apart by delta-s= distribCDFINV(1-PRC1)/NN1. Thus, the imprecision is about (MAX(1/NN1,PRC1). The lower the PRC1, and the higher the number of transit compartments in the expansion, the more precisely the convolution process is calculated. If you wish to minimize the computational expense, and still have a reasonable rough approximation to the distribution (and remember, the selected distribution function is in itself often arbitrary or empirically assessed and rarely associated with a mechanistic basis, so theoretical purity should not be a priority here), one could select NN1=10 and PRC1=0.01.
The user must have in his code the placement of the 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 section I.34 Probability Density Functions (NM742)), indices starting at position 2 should have the appropriate parameter values specified (the first index is for the variable of integration t, and is filled in by ddexpand, as needed):
; Define parameters to Weibull distribution
VG(2)=AA
VG(3)=BBThe above lines must be placed at beginning of $DES. The parameter locations for the distributions are ordered as follows, starting with index 2, using the parameter names listed in (Koch and Schropp 2015):
WEIBULLCDF: a,bGAMMACDF: a,bLOGNORMALCDF: b,sigmaGOMPMAKECDF(Gompertz-Makeham): a,b,cWEIBULLCCDF(combined Weibull): a1,b1,a2,b2,dUNIFORMCDF: a,b
When the above control stream is submitted to ddexpand (using the -dde option to the nmfe script), we get the resulting control stream, which carries out the afore-mentioned convolution of the distribution with the intended state variable, with the histogramming resolution of NN1, and up to the time determined by PRC1:
$DES
VG(2)=NU1
VG(3)=TAU1
APG_1_1=Y0
NNU_1=27
VG(1)=1.0-0.01
TTEND_1=WEIBULLCDFINV(VG)
KK_1=NNU_1/TTEND_1
DTT_1=1.0/KK_1
IF(TTEND_1==-1.0D+300) EXIT 2
VG(1)=DTT_1
"VG(11)=LOC(WEIBULLCDF)
VG(12)=-1.0D+300
DADT(2)=KK_1*(A(1)-A(2)-APG_1_1)
VG(13)=A(2)
DADT(3)=KK_1*(A(2)-A(3))
VG(14)=A(3)
...
DADT(28)=KK_1*(A(27)-A(28))
VG(39)=A(28)
VG(40)=-1.0D+300
SUM_1=TRANSITG_SDIM1(VG)
ADG_1_1=SUM_1+A(28)+APG_1_1
DADT(1)=KG*(1.0-ADG_1_1/YSS)*A(1)Sometimes you want to 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 APG_x_y. For example, in the following example, kin is to be convolved with a Weibull distribution (distrib_delay\smax_weibulldt.ctl):
;DDE
;NN1=100
;PRC=0.01
;DISTRIB=WEIBULLCDF(VG)
$PROBLEM LOGISTIC Distributed transit compartments on Smax Drug Responsive
Input function
…
$DES
VG(2)=AA
VG(3)=BB
c = A(1)/V
kin = kin0*(1.0+(Smax*c)/(SC50+c))
APG_KIN_1=KIN0
DADT(1) = -(Vmax*A(1))/(V50+A(1))
DADT(2) = kin – ADG_KIN_1
…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=100
;PRC=0.01
;DISTRIB=WEIBULL(VG)but generally results are better when using the CDF.
You may have a mixture of some states be general-distributed delayed (using APG_x_y and ADG_x_y), some states be discrete delayed (using AP_x_y and AD_x_y), and some states transit compartment delayed (using APT_x_y and ADT_x_y). The following is a portion of the epo problem, highlighting the delay modeling portion, in which ADG_4_1 is disributed delayed, and the others are discrete delayed.
;DDE
;NU1=10
;DISTRIB=WEIBULLCDF(VG)
$PROB EPO
$ABBR DERIV2=NO DERIV2=NOCOMMON DERIV1=NO
$INPUT ID AMT TIME DV EVID MDV CMT
$DATA test_epo.csv IGNORE=C
$SUBROUTINES ADVAN16 TOL=10 ATOL=10
$MODEL NCOMPARTMENTS=6
$PK
...
TP1 =THETA(17)
TP2 =THETA(18)
TRET =THETA(19)
TRBC =THETA(20)
NUG=5
; All Nuy are defined here, but only NU1 was needed and needs to be defined.
NU1=NUG
NU2=NUG
NU3=NUG
NU4=NUG
NU5=NUG
NU6=NUG
;TAUy
TAU1=TP1+TP2
TAU2=TP2
TAU3=TP1+TP2+TRET
TAU4=TP2+TRET
TAU5=TP1+TP2+TRET+TRBC
TAU6=TP2+TRET+TRBC
...
$DES
APG_4_1=RC0; Note the APG definition for distributed delay
AP_4_2=RC0
AP_4_3=RC0
AP_4_4=RC0
AP_4_5=RC0
AP_4_6=RC0
AP_5_1=RET0
AP_5_3=RET0
AP_5_5=RET0
AP_6_1=RBCM0
AP_6_3=RBCM0
AP_6_5=RBCM0
CC=A(1)/VP
AT=A(2)
RR=A(3)
RC=A(4)
RET=A(5)
RBCM=A(6)
; Note the accompanying ADG usage for distributed delay
X1=1+SMAX*ADG_4_1/(SC50+ADG_4_1)
X2=1+SMAX*AD_4_2/(SC50+AD_4_2)
X3=1+SMAX*AD_4_3/(SC50+AD_4_3)
X4=1+SMAX*AD_4_4/(SC50+AD_4_4)
X5=1+SMAX*AD_4_5/(SC50+AD_4_5)
X6=1+SMAX*AD_4_6/(SC50+AD_4_6)
X0=1+SMAX*RC0/(SC50+RC0)
I1=1-IMAX*(MCH*(AD_5_1+AD_6_1)-HB0)/(IC50+(MCH*(AD_5_1+AD_6_1)-HB0))
I3=1-IMAX*(MCH*(AD_5_3+AD_6_3)-HB0)/(IC50+(MCH*(AD_5_3+AD_6_3)-HB0))
I5=1-IMAX*(MCH*(AD_5_5+AD_6_5)-HB0)/(IC50+(MCH*(AD_5_5+AD_6_5)-HB0))
DADT(1)=KEPO-KON*CC*VP*RR+KOFF*RC*VP-(KEL+KPT)*CC*VP+KTP*AT
DADT(2)=KPT*CC*VP-KTP*AT
DADT(3)=KSYN-KON*CC*RR+KOFF*RC-KDEG*RR
DADT(4)=KON*CC*RR-(KOFF+KINT)*RC
DADT(5)=KIN*X1*X2*I1-KIN*X3*X4*I3
DADT(6)=KIN*X3*X4*I3-KIN*X5*X6*I5Although the discrete delays may be expanded by MOS (;MOS comment) or DDE (;DDE comment), expanding by ;DDE is generally more efficient.
13 Using Discrete Delay Solvers for Time Delays of any distribution, Using the ddexpand program (nm751) and Built-IN DDEFUNCG Routine
So far, transit compartments were used to approximate delay compartments for use in distributed delays. You can use the direct delay DDE algorithms to perform a distributed delay as well, which create perfect discrete delay evaluations of the past on state variables A(), instead of using the transit compartment method. The advantage to using the DDE solvers for distributional delays is that no additional differential equations are added to the user’s model system, and particularly with the built-in DDE solvers ADVAN16 to ADVAN18, are more efficient than the transit compartment distributed method. The user models direct delay solver method the same way you would distributed delays via transit compartments, 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)And of course you must use one of the ODE solvers that are capable of discrete delay solving (ADVAN9, ADVAN13, ADVAN14, ADVAN15, ADVAN16, ADVAN17, ADVAN18). Doing so will allow the use of the discrete DDE solvers on NN1 evenly spaced times t-i*delta_tau, for delay times i=1 to NN1, weighted according to the CDF of the distribution. 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, which would be NN1 in this example.
One down-side to using the built-in DDEFUNCG routine and non-dae discrete delay solvers is that variables such as KIN cannot be directly delayed, and require access via an equivalent state variable A(), in the manner we showed earlier. For example, you would add something like:
DADT(x)={expression for derivative of kin with respect to time}.In a manner that was shown earlier. In the example above, kin could be coded by the user as an additional state variable, A(3) (distrib_delay\smax_weibulldd.ctl):
;DDE
;NN1=100
;DISTRIB=WEIBULL(VG,2)
$PROBLEM LOGISTIC Distributed discrete delay on Smax Drug Responsive Input
function
; turn off second derivative assessments, sometimes even 1st derivatives if
only simulating
$ABBR DERIV2=NO DERIV2=NOCOMMON
$INPUT ID AMT TIME DV EVID MDV
$DATA smax_weibull.dat IGNORE=C
$SUBROUTINES ADVAN16 TOL=10 ATOL=10 OTHER=DDEFUNCG.F90
$MODEL NCOMPARTMENTS=3
$PK
CALLFL=-2
MXSTEP=2000000000
SMAX=THETA(1)
SC50=THETA(2)
KIN0=THETA(3)
AA=THETA(4)
BB=THETA(5)
Vmax = 2
V50 = 3
V=1
; Initial conditions
A_0(1)=0.0
A_0(2)=5.0
; KIN at time 0 actually contains the bolus dose contribution into A(1).
; So it is not kin0, but rather;
CCC = AMT/V
A_0(3)=kin0*(1.0+(Smax*ccc)/(SC50+ccc))
$DES
VG(2)=AA
VG(3)=BB
c = A(1)/V
; DERIVATIVE OF A WRT T, WHICH EQUALS V TIMES DERIVATIVE OF C WRT T
DADT1=-(Vmax*A(1))/(V50+A(1))
; DERIVATIVE OF KIN WRT C
DKINC=KIN0*SMAX*SC50/(SC50+C)/(SC50+C)
; The past, is when T<0, so before the bolus dose into A(1) is given, which
is:
APG_3_1=KIN0
DADT(1) = DADT1
; DERIVATIVE OF KIN WRT C TIMES DERIVATIVE OF C WRT T YIELDS DERIVATIVE OF
KIN WRT T
DADT(3) = DKINC*DADT1/V
; SUBMIT STATE VARIABLE A(3), WHICH REPRESENTS KIN, FOR WEIBULL-DISTRIBUTED
DELAY
DADT(2) = A(3) - ADG_3_1
$ERROR
CC = A(1)/V
kin = kin0*(1.0+(Smax*cc)/(SC50+cc))
A1=A(2)
KINCALC=A(3)
Y1=1.0
IPRED=A(2)
Y=IPRED*(1.0+EPS(1))
$THETA
0.5
10.0
0.8
7.0
6.5
$OMEGA (0.0 fixed)x5
$SIGMA (0.0 fixed)
$SIML (122345) ONLYSIM SUBP=1
;kin and kincalc will have identical values, as expected.
$TABLE TIME KIN KINCALC IPRED NOPRINT NOAPPEND FILE=smax_weibulldd.tab
ONEHEADER FORMAT=S1PE20.1313.1 Default quadrature behavior
By default, the ddexpand routine will establish the most suitable quadrature method for performing the convolution integration. The default is as follows:
If DISTRIB/=GAMMA and DISTRIB/=UNIFORM, then the Gaussian-Laguerre Gaussian Quadrature (GQ) method with shaping parameter SHP=-0.5 is used, in the gamma weight function X**SHP*exp(-x). This has been found to be the most robust and efficient for most distributions. The default PRC is 0.0001, unless it is specified with ;PRC.
If DISTRIB=UNIFORM, then the Gauss-Legendre GQ method is used, with lower bound and upper bound being that of the (a,b) parameters of the uniform distribution modeled.
If DISTRIB=GAMMA, then the Gaussian-Laguerre GQ method, with SHP=NU-1 (where NU is the shape parameter of the modeled Gamma distribution given by what the user specifies for VG(2)), if NU>=0.1. If NU<0.1, then this GQ method performs poorly, and instead the Gauss-Legendre method with domain range (0,1) is used on the quantile domain, using GAMMACDFINV for the re-mapping of the abscissas that the Gauss-Legendre algorithm produced. The default PRC is 0.0001, unless it is specified with ;PRC.
13.2 Overriding the default quadrature
Rarely does a user need to modify the default behavior. However, should it be necessary, the following Gaussian-Quadrature techniques are available, by specifying the algorithm as follows (to be placed after the ;DISTRIB line):
;ALGx=GQ({U/N/G},NU)where x refers to the xth index TAU of the total model. By default, x is the most recent index specified by ;NNx. GQ refers to using a Gaussian Quadrature method. The arguments can be given in any order, as there is no risk of confusion. The Weighting function of the GQ method could be:
U=uniform=Gauss-Legendre
N=normal=Gauss-Hermite
G=Gamma=Gauss-LagerreThe second argument NU is optional, and would be valid only when G weighting function is used. In this case:
If NU is not specified, the default alpha is determined in the manner described earlier.
If NU is specified >0, then it is used.
If NU>=1000, then the model’s NU=VG(2) is used. Use NU>=1000 only if DISTRIB=GAMMA.
If NU=0, then when DISTRIB/=GAMMA, the default SHP=-0.5 (NU=0.5) is used. When DISTRIB=GAMMA, then the model parameters are used, where SHP=VG(2)-1. In this case, the PRC value is not used for specifying integral range, and instead the scaling using SCL=VG(3) is done on the default time positions (and therefore integral range) of the GQ method. The relationship between VG(3) and the resulting effective PRC is complex, and may result in insufficient integral range and accuracy.
If U is selected, then the Gauss-Legendre->DISTRIBCDFINV(abscissas) process is used, where abscissas is the vector of positions calculated by Gauss-Legendre.
If N is selected, then the Gauss-Hermite->DISTRIBCDFINV(exp(abscissas)) is used. This provides poor results, even for the lognormal distribution, so it should not be used.
From experience, it appears that the default GQ method (so, not specifying ALG), with NNx=100, tends to very often provide a precision of the convolution equivalent to RTOL=0.01, ATOL=0.01, that is, 1% imprecision. An exception is for GAMMA with NU<=0.1, in which case NNx=150 may be required to get a 1% imprecision. A few trial and error simulations with varying NNx is sufficient to verify the imprecision, and that NNx can then be used for the entire analysis.
13.3 Newton-Cotes quadratures
The Newton-Cotes quadratures may also be used, but they are not as efficient as the GQ methods. You may want to use Newton-Cotes to have more standard histogrammic representation of the convolution of limited resolution (NNx=10, say), to purposely have a discretized rendition of the convolution, perhaps to replicate a result from the literature, in which the convolution analysis may have been simplified (as mentioned earlier, if there is no good mechanistic reason to have a specific distribution, but merely to represent some empirical average life-span with some degree of variability, this could be sufficient). This syntax is as follows:
ALGx=NC(n,{T/Q},PRC,{P/C})The arguments may be given in any order, as there is no risk of confusion. The Newton Cotes algorithms are as follows:
n=1: rectangular steps of integration
n=2: trapezoidal rule
n=3: Simpson's rule
n=4: Simpson's 3/8 rule
n=5: Bode's rule
6<=n<=35: Newton-Cotes{T/Q}: Select T for equally spaced time values throughout the integration, select Q for equally spaced quantile values (that is, for which CDF values are equally spaced (find taui such that CDF(taui)=i/m)). Default is Q.
PRC: Value between 0 and <1, to indicate that the integration should be performed from 0 to the 1-PRC quantile position, similar to what is done for the transit compartment distributed delay method described in the previous section. Default is 0.1/NNx for Newton-Cotes. The PRC may also be specified on a separate line with ;PRC.
{P/C}: Select to use CDF (C) of the distribution, or PDF (P). Default is C, and generally gives the best results. P can only be selected with T. Selecting P is risky, for pdf’s that have a non-defined value, such as NU<1 for a GAMMA density.
When Q is selected and n is not specified, then default is n=5, as it appears to be the most accurate for most problems. When T is selected and n is not specified, then default n=3.
So:
ALGx=NCchooses n=5, Q, C.
ALGx=NC(T,0.001)chooses 3, T, C, and PRC=0.001.
ALGx=NC(2,T)selects trapezoidal rule on equally spaced times, using CDF in the integrand.
ALGx=NC(3,T,P)Selects Simpson’s rule on equally spaced times, using pdf in the integrand.
13.4 Equilibrium compartment method with dae DDE solver
We can alternatively use the equilibrium compartment method (using a dae dde solver) that we had shown earlier for a simple discrete delay problem, but now apply it for distributed delay as follows (distrib_delay\smax_weibulldd_dae.ctl):
;DDE
;NN1=100
;DISTRIB=WEIBULLCDF(VG,2)
$PROBLEM LOGISTIC Distributed discrete delay on Smax Drug Responsive Input
function using dae
$ABBR DERIV2=NO DERIV2=NOCOMMON
$INPUT ID AMT TIME DV EVID MDV
$DATA smax_weibull.dat IGNORE=C
$SUBROUTINES ADVAN17 TOL=10 ATOL=10
$MODEL NCOMPARTMENTS=3
COMP(COMP1)
COMP(COMP2)
COMP(COMP3 EQUILIBRIUM)
$PK
CALLFL=-2
MXSTEP=2000000000
SMAX=THETA(1)
SC50=THETA(2)
KIN0=THETA(3)
AA=THETA(4)
BB=THETA(5)
Vmax = 2
V50 = 3
V=1
; Initial conditions
A_0(1)=0.0
A_0(2)=5.0
; KIN at time 0 actually contains the bolus dose contribution into A(1).
; So it is not kin0, but rather;
$DES
VG(2)=AA
VG(3)=BB
APG_3_1=kin0
DADT(1) = -(Vmax*A(1))/(V50+A(1))
DADT(2) = A(3) - ADG_3_1
$AESINIT
INIT=0
$AES
c = A(1)/V
kin = kin0*(1.0+(Smax*c)/(SC50+c))
E(3)=A(3)-kin
$ERROR
A1=A(2)
Y1=1.0
IPRED=A(2)
Y=IPRED*(1.0+EPS(1))
$THETA
0.5
10.0
0.8
7.0
6.5
$OMEGA (0.0 fixed)x4
$SIGMA (0.0 fixed)
$SIML (122345) ONLYSIM SUBP=1
$TABLE TIME KIN IPRED NOPRINT NOAPPEND FILE=smax_weibulldd_dae.tab ONEHEADER
FORMAT=S1PE20.1313.5 Adding an additional discrete time shift to the delayed distribution (nm760)
For a distributed delay pertaining to any APG_x_y, you can also add a discrete time shift TSHIFT_x, or discrete delay, to the distributed delay equation. For example, in the following code:
;DDE
;NN1=10
;PRC1=0.01
;DISTRIB1=GAMMA(VG,2)
…
$PK
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
MU_4=THETA(4)
MU_5=THETA(5)
MU_6=THETA(6)
KG=EXP(MU_1+ETA(1))
Y0=EXP(MU_2+ETA(2))
YSS=EXP(MU_3+ETA(3))
KK=EXP(MU_4+ETA(4))
NU1=EXP(MU_5+ETA(5))
TSHIFT=EXP(MU_6+ETA(6))
A_0(1)=Y0
VG(2)=NU1
VG(3)=KK
$DES
VG(2)=NU1
VG(3)=KK
TSHIFT_1=TSHIFT
APG_1_1=Y0
DADT(1)=KG*(1.0-ADG_1_1/YSS)*A(1)
…The distributed delay is effectively TSHIFT+GAMMA(T…). If TSHIFT is 10 hours, then the delay ADG_1_1 is state A(1) delayed by 10 hours, followed by distribution of the state variable in accordance with the Gamma distribution. This is equivalent to using a delay distribution of function Gamma(T-TSHIFT…), so its domain is from T+TSHIFT to T=infinity.
14 General Quadrature Integration Using ddexpand, and Routines GQ and NC
In the previous section, we used a packaged routine called DDEFUNCG that efficiently evaluates quadrature integration of distributed time delays using pre-packaged 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 used. This allows for more rapid evaluation. Furthermore, the coding of the convolution process is done for you. The down-side is that 1) only convolution integrations using pre-packaged distributions could be used as the kernel function, 2), the past (AP_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 the afore-mentioned limitations, is to have the user code within NMTRAN the specific equations to be integrated, and using the DOWHILE feature, along with calling the quadrature routines GQ or NC directly.
Revisiting the smax_weibulldd 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, ..\examples\distrib_delay\smax_weic.ctl, but requiring the user to code more of the equations:
;DDE
$PROBLEM Discrete Delay on Smax Drug Responsive Input function
; turn off second derivative assessments, sometimes even 1st derivatives if
only simulating
$ABBR DERIV2=NO DERIV2=NOCOMMON DERIV1=NO
$INPUT ID AMT TIME DV EVID MDV
$DATA smax_gammat.dat IGNORE=C
$SUBROUTINES ADVAN16 TOL=14 ATOL=14
$MODEL NCOMPARTMENTS=2
$PK
CALLFL=-2
MXSTEP=2000000000
SMAX=THETA(1)
SC50=THETA(2)
KIN0=THETA(3)
AA=THETA(4)
BB=THETA(5)
Vmax = 2
V50 = 3
V=1
TEND=BB*(-LOG(0.0001))**(1.0D+00/AA) ; Have TEND cover most of the integral
; Initial conditions
A_0(1)=0.0
A_0(2)=0.0
$DES
DADT(1)=0.0
DADT(2)=0.0
SUM=0.0
SUMW=0.0
;DOC 25 TS CS GQ(G,0.5,TEND)
TAU1=TS
AP_1_1=0.0
LL=AA/BB*((TS/BB)**(AA-1.0))*EXP(-((TS/BB)**AA))
CDELAY=AD_1_1/V
FF=KIN0*(1.0+(SMAX*CDELAY)/(SC50+CDELAY))
SUM=SUM+CS*LL*FF
SUMW=SUMW+CS*LL
;ENDDOC
ADC=SUM/SUMW
C=A(1)/V
DADT(1)=-(VMAX*A(1))/(V50+A(1))
KINW=KIN0*(1.0+(SMAX*C)/(SC50+C))
DADT(2)=KINW-ADC
$ERROR
CC=A(1)/V
A1=A(2)
Y1=1.0
IPRED=A(2)
Y=IPRED*(1.0+EPS(1))
kin = kin0*(1.0+(Smax*cc)/(SC50+cc))
$THETA
0.5
10.0
190.0
7.0
6.5
$OMEGA (0.0 fixed)x4
$SIGMA (0.0 fixed)
$SIML (122345) ONLYSIM SUBP=1
$TABLE TIME KIN IPRED NOPRINT NOAPPEND FILE=smax_weic.tab ONEHEADER
FORMAT=S1PE20.13The ;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 next item is the number of points for the integration, which is 25 in this example. The next item is the variable of quadrature integration, which is TS in this example. The next item is the coefficient of quadrature integration, which is CS in this example. The next item is the quadrature method to be used, along with arguments. GQ or NC may be selected, with the following arguments:
GQ(L,A,B):
L=U=gaussian-legendre (uniform weighting) quadratureA=beginning Time,B=ending timeL=G=gaussian-laguerre (gamma weigthing) quadratureA=NU,B=ending timeTEND. Beginning time is always 0.L=N=gaussian-Hermite (normal function, times then exponentiated).A=beginning Time,B=ending time
NC(x,A,B):
X=x-point Newton-Cotes methodA=beginning Time,B=ending time
Notice how the variable of integration TS and coefficient of integration CS are used in the accumulator variables SUM and SUMW. The FF variable is kin calculated using the state variable A(1) delayed by an amount 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.
If you are specifically convolving with the GAMMA distribution, you can use the Gauss-Laguerre’s gamma function weighting to its full advantage as shown in the following example (..\examples\distrib_delay\logits_ddg.ctl):
;DDE
$PROBLEM LOGISTIC
$ABBR DERIV2=NO DERIV2=NOCOMMON ; DERIV1=NO
$INPUT ID AMT TIME PRDV DV EVID MDV
$DATA logistic6.csv IGNORE=C
$SUBROUTINES ADVAN16 TOL=6 ATOL=6
$MODEL NCOMPARTMENTS=1
$PK
MXSTEP=2000000000
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
MU_4=THETA(4)
MU_5=THETA(5)
KG=EXP(MU_1+ETA(1))
Y0=EXP(MU_2+ETA(2))
YSS=EXP(MU_3+ETA(3))
KK=EXP(MU_4+ETA(4))
NU1=EXP(MU_5+ETA(5))
A_0(1)=Y0
$DES
ADC=0.0
;DOC1 10 TS CS GQ(GAM,NU1,KK)
TAU1=TS
AP_1_1=Y0
ADC=ADC+CS*AD_1_1
;ENDDOC1
DADT(1)=KG*(1.0-ADC/YSS)*A(1)
$ERROR
A1=A(1)
Y1=1.0
IPRED=A(1)
Y=IPRED*(1.0+EPS(1))
$THETA
-1.61764E+00
-2.40835E-02
2.29145E+00
-0.5
1.098612289
$OMEGA BLOCK(3) 6.86871E-03
8.45335E-03
6.59591E-04 7.39132E-03
1.01727E-03 -1.73680E-03
$OMEGA 0.01 0.01
$SIGMA
3.23154E-03
;$SIML (122345) SUBP=1
;$EST METHOD=IMP AUTO=1 NOABORT PRINT=1
$EST METHOD=CHAIN FILE=temp14.ext ISAMPLE=-1000000000 TBLN=1 NSAMPLE=0
$EST METHOD=ITS INTERACTION NOHABORT SIGL=5 MCETA=10 NSIG=2 PRINT=1 NITER=0
CTYPE=3 OPTMAP=0 ETADER=0 FNLETA=0
$EST METHOD=ITS INTERACTION NOHABORT SIGL=5 MCETA=10 NSIG=2 PRINT=1 NITER=0
CTYPE=3 OPTMAP=1 ETADER=2 FNLETA=0Make sure you specify “GAM” instead of “G” for the GC function. Notice that the CS coefficients already contain the gamma distribution weighting, appropriate specifically for parameters NU1 and KK, so that the density need not be multiplied into the summand. Note also that the appropriate TEND will be determined based on NU1 and KK.
Suppose it is desired to create quadrature positions along quantile positions of the distribution. This might be particularly useful if there may be singularities in the the pdf. The following does this for the gamma density (distrib_delay\logits_ddcdfinv.ctl):
;DDE
$PROBLEM LOGISTIC
$ABBR DERIV2=NO DERIV2=NOCOMMON ; DERIV1=NO
$ABBR FUNCTION GAMMACDFINV(VG,10)
$INPUT ID AMT TIME PRDV DV EVID MDV
$DATA logistic6.csv IGNORE=C
$SUBROUTINES ADVAN16 TOL=6 ATOL=6
$MODEL NCOMPARTMENTS=1
$PK
MXSTEP=2000000000
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
MU_4=THETA(4)
MU_5=THETA(5)
KG=EXP(MU_1+ETA(1))
Y0=EXP(MU_2+ETA(2))
YSS=EXP(MU_3+ETA(3))
KK=EXP(MU_4+ETA(4))
NU1=EXP(MU_5+ETA(5))
DEN=KK/EXP(GAMLN(NU1))
A_0(1)=Y0
$DES
VG(2)=NU1
VG(3)=KK
ADC=0.0
;DOC1 10 TSS CS GQ(U,0.0,1.0)
VG(10)=II_1
VG(1)=TSS
TS=GAMMACDFINV(VG)
TAU1=TS
AP_1_1=Y0
ADC=ADC+CS*AD_1_1
;ENDDOC1
DADT(1)=KG*(1.0-ADC/YSS)*A(1)
$ERROR
A1=A(1)
Y1=1.0
IPRED=A(1)
Y=IPRED*(1.0+EPS(1))
$THETA
-1.61764E+00
-2.40835E-02
2.29145E+00
-0.5
1.098612289
$OMEGA BLOCK(3) 6.86871E-03
8.45335E-03
6.59591E-04 7.39132E-03
1.01727E-03 -1.73680E-03
$OMEGA 0.01 0.01
$SIGMA
3.23154E-03
$EST METHOD=CHAIN FILE=temp14.ext ISAMPLE=-1000000000 TBLN=1 NSAMPLE=0
$EST METHOD=ITS INTERACTION NOHABORT SIGL=5 MCETA=10 NSIG=2 PRINT=1 NITER=0
CTYPE=3 OPTMAP=0 ETADER=0 FNLETA=0
$EST METHOD=ITS INTERACTION NOHABORT SIGL=5 MCETA=10 NSIG=2 PRINT=1 NITER=0
CTYPE=3 OPTMAP=1 ETADER=2 FNLETA=0In the above example, the Gaussian-Legendre guadrature is selected, with returned time values TSS scaled between 0.0 and 1.0, the domain of a CDF inverse (or range of the CDF). These TSS are transformed into time positions TS using the GAMMACDFINV function. Because GAMMACDFINV is an expensive calculation, the routine has the option to store up to 1000 values, so every unique value needs to be calculated only once, and then recalled upon repeated use. Each unique result needs to be given a unique index, and a convenient one is the DOWHILE index II_y (ddexpand creates II_y for a given ;DOCy). IF NDIM of the VG() vector is at least 10, then the tenth argument (VG(10)) must contain an index number. Under these conditions, the GAMMACDFINV routine will test if any of the input variables for a given index II_y changed, and will re-calculate if it has changed, or recall from the internal array if the input variables have not changed. Other CDFINV functions such LOGNORMALCDFINV, GOMPMAKECDFINV, and WEIBULLCCDFINV, also have this efficiency method, with the unique index expected to be entered in the 10th element of the input array.
The integral need not be a classic convolution with a distribution. The following example shows the considerably versatility of this method. The problem ..\examples\distrib_delay\lympho.ctl is as follows:
;DDE
$SIZES MAXNRDS=2 PAST_SIZE=20000
$PROBLEM DDDE Model of Lymphocyte Trafficking
; turn off second derivative assessments, sometimes even 1st derivatives if
only simulating
$ABBR DERIV2=NO DERIV2=NOCOMMON DERIV1=NO
$ABBR FUNCTION GAMMACDF(GCDF,4)
$INPUT ID AMT TIME DV EVID MDV
$DATA lympho.dat IGNORE=C
$SUBROUTINES ADVAN16 TOL=6 ATOL=6
$MODEL NCOMPARTMENTS=3
$PK
CALLFL=-2
MXSTEP=2000000000
AB0=THETA(1)
SHAZ0=THETA(2)
KOUT=THETA(3)
U1SHAZ0=THETA(4)
NU=THETA(5)
IMAX=THETA(6)
IC50=THETA(7)
CL=THETA(8)
VP=THETA(9)
U2=(1-U1SHAZ0)/SHAZ0
KIN=AB0*KOUT*U2*SHAZ0
GNU=EXP(GAMLN(NU))
QQ=(U1SHAZ0)**(1.0/NU)
BETA=U2*QQ/(1.0-QQ)
; Initial conditions
A_0(1)=0.0
A_0(2)=AB0
A_0(3)=KIN/U2
TEND=-LOG(0.0001)/U2 ; TEND is calculated to be when exp(-u2*t)<=0.001, so
should be enough
$DES
DADT(1)=0.0
DADT(2)=0.0
DADT(3)=0.0
CP=A(1)/VP
EMAX=1.0-IMAX*CP/(IC50+CP)
ADC1=0.0
ADC2=0.0
;DOC1 30 TA CA GQ(U,0.0,TEND)
IUU=0.0 ; IUU=INTEGRATED U=U1+U2
DEN=GNU
;DOC2 30 TZ CZ GQ(U,0.0,TA)
TAU2=TA-TZ
IF(TAU2<0.0) TAU2=0.0
AP_1_2=0.0
CPD2=AD_1_2/VP
EMAX2=1.0-IMAX*CPD2/(IC50+CPD2)
U1TZ=BETA*((BETA*TZ)**(NU-1.0))*EXP(-BETA*TZ)
DEN=DEN-U1TZ
GCDF(1)=TZ
GCDF(2)=NU
GCDF(3)=BETA
DENZ=GNU-GAMMACDF(GCDF)
IUU=IUU+CZ*(EMAX2*U1TZ/DENZ+U2)
;ENDDOC2
GCDF(1)=TA
GCDF(2)=NU
GCDF(3)=BETA
DENZ2=GNU-GAMMACDF(GCDF)
U1=BETA*((BETA*TA)**(NU-1.0))*EXP(-BETA*TA)/DENZ2
SHAZ=EXP(-IUU)
TAU1=TA
AP_2_1=AB0
CPD=AD_2_1
ADC1=ADC1+CA*CPD*U1*SHAZ
ADC2=ADC2+CA*CPD*(U1*EMAX+U2)*SHAZ
;ENDDOC1
DADT(1)=-CL*CP
DADT(2)=KIN-KOUT*A(2)+KOUT*EMAX*ADC1
DADT(3)=KOUT*A(2)-KOUT*ADC2
$ERROR
CC=A(1)/VP
Y1=1.0
IPRED=CC
A2=A(2)
A3=A(3)
Y=IPRED*(1.0+EPS(1))
$THETA
31.4
48.0
0.137
0.5
1.0
1.0
4.34
9471.0
70961.0
$OMEGA (0.0 fixed)x5
$SIGMA (0.0 fixed)
$SIML (122345) ONLYSIM SUBP=1
$TABLE TIME IPRED A2 A3 NOPRINT NOAPPEND FILE=lympho.tab ONEHEADER
FORMAT=S1PE20.13It is beyond the present scope to explain the above model, but there are several features used here that have general applicability. Notice that Do while loop 2 (;DOC2) is nested within DO while loop 1 (;DOC1). When ddexpand expands these comment directives, it defines index variable II_y for DO While loop DOCy, and termination index NN_y. To distinguish the various do while activities, always have a unique y value for every ;DOCy definition in the control stream. Notice that each ;DOC has its own integration variable and coefficient of integration variable.
One caveat should be noted when using ;DOC (or more generally, DOWHILE). Only accumulator variables (those defined as value=value+…) have their analytical derivative dependencies carried over properly from one DOWHILE iteration to the next. So any dependency on a value of the previous iteration, such as TSPRE AND CSPRE in this example:
;DOC1 13 TS CSS NC(3,0.0,TEND)
…
DELT=TS-TSPRE
DELC=CS-CSPRE
IF(DELT/=0.0) ADC=ADC+CSS*DELC*AD_1_1/DELT
TSPRE=TS
CSPRE=CS
;ENDDOC1will not have their derivatives properly calculated. Model instead as follows:
TSPRE=0.0
CSPRE=0.0
;DOC1 13 TS CSS NC(3,0.0,TEND)
…
DELT=TS-TSPRE
DELC=CS-CSPRE
IF(DELT/=0.0) ADC=ADC+CSS*DELC*AD_1_1/DELT
TSPRE=DELT+TSPRE
CSPRE=CSPRE+DELC
;ENDDOC115 summary: A utility to provide summary statistics of a series of table outputs (NM751)
After creating a series of table outputs generated by a $SIML statement, it may be desired to obtain summary statistics of the table items. The summary utility program can provide these summary statistics. For example, suppose you generate the table outputs from a control stream as follows:
$SIM (4442223) SUBPROB=1000
…
$TABLE ID TIME … IPRED … NOPRINT NOAPPEND FILE=table.tabYou may want to obtain quantile values of these 1000 generated samples of each table item. The syntax of executing this program is as follows (the third argument, a FORTRAN style format, is optional):
summary table.tab summary.tab "1X,1PG15.8"In the above example, the resulting summary.tab file will contain quantile position values for each of the table items, for example (RSTD=relative standard deviation):
TIME: Mean STD RSTD 2.50% 97.50%
Row 0.0000 0.0000 0.0000 0.0000 0.0000
36.563 6.32440E-02
1 0.13934 5.09463E-02 24.814 4.0921 0.26127
2 7.1222 1.7673 24.814 4.0921 11.134
3 7.1222 1.7673 24.814 4.0921 11.134
4 7.1222 1.7673 24.177 82.548 11.134
5 155.81 37.671 199.97
6
IPRED: Mean STD RSTD 2.50% 97.50%
Row 0.0000 0.0000 0.0000 0.0000 0.0000
1.0839 12.732 0.84534 1.3868
1 7.9587 0.13800 26.708 4.5215 12.773
2 7.9587 2.1256 26.708 4.5215 12.773
3 7.9587 2.1256 26.708 4.5215 12.773
4 2.1256 39.702 0.37273 1.1902
5 0.52453
6 0.2082516 lpcc: LinePrinter Controls Converter Program (NM75)
The utility program lpcc.f90 in the ..\util directory will convert the line printer control characters in the NONMEM report file into their proper actions. The line printer control characters are in the first character in each line, and at one time was interpreted by old 132 column line-printers.
The conversion is acted upon only between the lines
NM-TRAN MESSAGESand just before
Stop Time:as the lines before and after these do not use column 1 as a line-printer control column.
The code command characters are interpreted by lpcc as follows:
| Code Option | Acting on Control character | Action |
|---|---|---|
| A0 | + |
Advance 0 lines (overprint on (merge with) previous line) |
| A1 | space | Advance 1 line (no action is taken) |
| A2 | 0 |
Advance 2 lines (insert a blank line) |
| A3 | - |
Advance 3 lines (insert 2 blank lines). Not used in NONMEM |
| AS | Remove action character column once code is acted upon (But will not remove character in col 1 if not one of the above) | |
| FF | 1 |
Advance page (insert Form feed character) |
| FFA | 1 |
Advance page and add blank line (insert Form feed with a next line) |
| FFA2 | 1 |
Advance page and add 2 blank lines |
| FF1 | 1 |
Replace '1' as ' ' (interpret FF as A1) |
| FF2 | 1 |
Replace '1' as '0' (interpret FF as A2) |
| FF3 | 1 |
Replace '1' as '-' (interpret FF as A3) |
| ALL | FF, A0, A1, A2, A3 | |
| ALLS | FF, A0, A1, A2, A3, AS |
After an action, the control character is replaced with a space. That space will be printed unless AS is requested. So AS shifts the line one column to the left.
An N followed by the code prevents that action from occurring.
Examples (options are acted upon in order):
lpcc example1.res example1.lst FF A0Only replace '1' in column 1 with form feed, and merge line with '+' with previous line.
lpcc example1.res example1.lst ALL FFAAct on all codes, and replace '1' with FF/next line.
lpcc example1.res example1.lst ALLS FFAAct on all codes, replace '1' with FF/next line, and shift line one column to the left.
lpcc example1.res example1.lst ALLS NFFAct on all codes except '1' (FF), and shift line one column to the left.
17 lblr: Insert User-Defined Labels in Additional Files (NM75)
The label substitution process (NOSUB=0 (Default) (NM74) in section I.32 General New Options for $ESTIMATION Record (NM73), and section I.8 Expansions on Abbreviated and Verbatim Code and Other Items) is normally done only in the NONMEM report file, but not in the additional files, such as the .ext, .cov, .coi, .phi, etc. These additional files need often be read by 3rd party software, which may rely on traditional names for THETAs, OMEGAs, SIGMAs, phi, eta, etc.
Nonetheless, one may wish to make label substitutions on these additional files, which can be done after the NONMEM analysis, using the label-replace utility, lblr.
The lblr utility is used as follows:
lblr myresult.ext myresult.res myresult_new.ext {S/T/C} {0/1/2}where
myresult.ext= input file to have label substitutions done.myresult.res= NONMEM report file containing label substitution patterns, or any file containing label substitution patterns such as:
THETA(1)=THETA(CL)
ETA(2)=ETA(ETV1)
EPS(1)=EPS(RSW)myresult_new.extis the new file with substitutions.
17.1 Option 1 settings
S=Expand and adjust with spaces (default)T=Do not expand, and truncate label if necessaryC=Expand, and compress so there are no spaces (suitable only for comma delimited files)
17.2 Option 2 settings
0=Full name of parameter type (THETA/OMEGA/SIGMA/ETA/EPS/PHI/ETC/PHC) (default)1=two letter truncation of parameter type name (TH/OM/SG/ET/EP/PH/EC/PC)2=single letter truncation of parameter type name (T/O/S/E/P/H/C/C)
Truncated parameter type name (1 or 2) works well with T, when no label expansion is permitted, and you want to make the label as compact as possible.
If the user specified NOSUB=1, then the THETAs, OMEGAs/etas, and SIGMAs/eps will not have been label substituted in the report file during the NONMEM execution. The user may wish to make the label substitution afterword with lblr. For example,
lblr myresult.res myresult.res myresult_new.res T 2Notice in this example that the report file myresult.res serves as the input file, and the file from which to get label substitution pattern. If you also want to perform an lpcc, do lblr first, then the lpcc utility.
18 csvalign: Pad a csv File with Spaces to align columns (NM75)
The utility program csvalign.f90 in the ..\util directory will pad a csv file with spaces so that the columns align, and more easily readable.
csvalign source.csv output.csv {number of header lines to exclude}For example, consider file 501.csv:
C,DataDesc:Problem3,,,,,,,
C,ID,TIME,DV,AMT,RATE,WT,AGE,SEX
.,1,0,0,1000,200,58.4,51,1
.,1,1,6.47,0,0,58.4,51,1You may not want to include the header line in considering the maximum spacing needed, so exclude it:
csvalign 501.csv 504_new.csv 1The result is:
C,DataDesc:Problem3,,,,,,,
C ,ID ,TIME ,DV ,AMT ,RATE ,WT ,AGE ,SEX
. ,1 ,0 ,0 ,1000 ,200 ,58.4 ,51 ,1
. ,1 ,1 ,6.47 ,0 ,0 ,58.4 ,51 ,119 c_add: Add Descriptor line and comment (C) first column to a data file (NM760)
The utility program c_add.f90 in the ..\util directory will add a descriptor line and command (C) first column to a data file, compatible with what PDx-Pop expects, and which can be read with a
$DATA … IGNORE=COr
$DATA … IGNORE=@NMTRAN control stream record.
19.1 Usage
c_add input.csv output.csv s/,/tDelimiter can be comma (,), t for tab, or s for space. Optional if delimiter is space.
19.2 Example
Input file:
ID,TIME,AMT,WT,APGR,DV,EVID,MDV
1,0,25,1.4,7,.,1,1
1,2,.,1.4,7,17.3,0,0
1,12.5,3.5,1.4,7,.,1,1Output file:
C,Data Desc: Place data description here
C,ID,TIME,AMT,WT,APGR,DV,EVID,MDV
.,1,0,25,1.4,7,.,1,1
.,1,2,.,1.4,7,17.3,0,0
.,1,12.5,3.5,1.4,7,.,1,120 priorget: Transfer Results of an Analysis to NMTRAN Prior Information (NM75)
This utility extracts results from a raw output file (.ext) and .cov file into $PRIOR NWPRI code, for easy transfer of results of a previous analysis into informative prior values for a subsequent analysis.
20.1 Usage
priorget root.ext root.cov myprior.ctl PRIORTYPE TBLN s/,/twhere root.ext and root.cov are input the files from a previous analysis.
PRIORTYPE=
Nif the previous analysis (source of prior information) was maximum likelihood/nlme (FOCE/Laplace/IMP/SAEM/ITS)Bif the previous analysis was Bayesian (BAYES/NUTS)Sif the previous analysis was Bayesian, and you want theSCALE()option of$OMEGAP/$SIGMAPused instead of modifying the matrix elements (letting NMTRAN scale the elements).Iif the previous analysis was Bayesian, and you want to use the inverse of the average of the inverse ofOMEGA/SIGMAof the previous analysis, and their variances for degrees of freedom assessment. According to statistical theory, the inverse of OMEGAS/SIGMAS from a Bayesian analysis are Wishart distributed, and they tend to be more similar in value to the mode of OMEGAs estimated in nlme analyses. See section I.43 A Note on Setting up Prior Information.
TBLN is the table number, as these input files may have several tables to choose from.
Delimiter can be ,, t for tab, or s for space. Optional if delimiter is space.
The order of the parameters in the input files must be in the default TSOL (lower triangular) order.
20.2 Example
priorget example3.ext example3.cov prior3.ctl N 5 sextracts information from table 5 of the population mixture problem example3, and places it in NMTRAN $PRIOR control stream format in the output file prior3.ctl:
$PRIOR NWPRI
$THETAP
( 4.24157E+00 FIXED)
(-2.29640E+00 FIXED)
( 4.26189E+00 FIXED)
(-6.75464E-01 FIXED)
( 6.66665E-01 FIXED)
$SIGMAP BLOCK(1)
1.02333E-02 FIXED
$OMEGAP BLOCK(2)
4.61841E-02 FIXED
-8.57594E-03 2.72391E-02
$OMEGAP BLOCK(2)
4.93083E-02 FIXED
-1.10257E-02 5.99622E-02
$THETAPV BLOCK(5) 6.00993E-04 7.40750E-04
2.40775E-04 FIXED 2.38694E-08
-5.00646E-05 1.51280E-04
1.10075E-07 1.90557E-08 5.12763E-04
-2.07126E-09 1.87141E-08 -1.13186E-04
-1.01416E-09 5.36580E-09 2.29895E-09
;DF from OMEGA(1,1): 184.063854545076
;DF from OMEGA(2,2): 162.626638917147
$OMEGAPD 162.626638917147 FIXED
;DF from OMEGA(3,3): 92.4955765039376
;DF from OMEGA(4,4): 99.3586543336461
$OMEGAPD 92.4955765039376 FIXED
$SIGMAPD 1780.99669390981 FIXEDYou can then patch (or “include”) this code into your control stream for the subsequent analysis. Any fixed THETAs will have $THETAPV variances with value 0, so you will need to modify this as needed. The degrees of freedom for OMEGAs and SIGMAs ($OMEGAPD, $SIGMAPD) are calculated according to the formula described in section I.43 A Note on Setting up Prior Information. Be aware that the degrees of freedom may fall below the block dimension, or be above the total number of subjects, and the user may wish to modify the results accordingly before using them in an analysis. Results from priorget should always be treated as informational first, and then inspected and modified according to your intended purpose.
21 nmtemplate Utility Program (NM73)
The utility program nmtemplate in the ..\util directory will perform variable substitution on appropriately tagged control stream template files, and produce executable control stream files. The syntax is as follows:
nmtemplate source-template-file destination-file var1=val1 var2=val2 var3=val3 …where var1=val1 is the variable name, and value to substitute in the template file. The variable var1 must in turn appear as <var1> in the template file, and is case sensitive. For example, consider the template file ..\util\nmtemp.nmt:
$PROB RUN# Example 1 (from samp5l)
$INPUT C SET ID JID TIME DV=CONC AMT=DOSE RATE EVID MDV CMT CLX V1X QX V2X SDIX SDSX
$DATA nmtemp2.csv IGNORE=C ACCEPT=(ID.EQ.<NMID>)
$SUBROUTINES ADVAN3 TRANS4
$PK
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
MU_4=THETA(4)
CL=DEXP(MU_1+ETA(1))
V1=DEXP(MU_2+ETA(2))
Q=DEXP(MU_3+ETA(3))
V2=DEXP(MU_4+ETA(4))
S1=V1
$ERROR
IPRED=F
Y = F + F*EPS(1)
; Initial values of THETA
$THETA <TH1> <TH2> <TH3> <TH4>
$OMEGA BLOCK(4)
0.15
0.01 0.15
0.01 0.01 0.15
0.01 0.01 0.01 0.15
$SIGMA
(0.06 )
$ETAS (0)x4
$EST METHOD=1 INTERACTION FNLETA=2 MAXEVAL=0
$TABLE ID TIME DV IPRED CMT EVID MDV ETA1 ETA2 ETA3 ETA4 NOAPPEND NOPRINT NOTITLE FILE=nmtemp.tabNote that <NMID> is to be replaced with a particular NONMEM ID number by nmtemplate, and the <THX> are to be replaced with specific values of THETAs:
nmtemplate nmtemp.nmt nmtemp.ctl NMID=47 TH1=1.7 TH2=1.4 TH3=0.8 TH4=2.0The resulting file nmtemp.ctl will have the various values substituted into the various <> placeholders, and is ready to be read by NMTRAN:
nmfe75 nmtemp.ctl nmtemp.resIn the above nmtemp.nmt example, because FNLETA=2, then NONMEM will simply evaluate the IPRED values using the inputted etas from the $ETAS record without performing an estimation. Another example template file is example6.nmt listed in the ..\util directory, that you may inspect for other ideas.
Actually, nmtemplate is a general variable substitution program, and can process any text file in the manner shown above. Consider a FINEDATA control stream file template (..\util\nmtemp.fnt):
$INPUT C SET ID JID TIME DV=CONC AMT=DOSE RATE EVID MDV CMT CLX V1X QX V2X
SDIX SDSX
$DATA nmtemp.csv IGNORE=C
$FINEDATA AXIS=TIME(LIN) TSTOP=<TSTOP> TSTART=<TSTART> NEVAL=<NEVAL>
FILE=nmtemp2.csvin which the tstart, tstop, and neval parameters are to be inserted:
nmtemplate nmtemp.fnt nmtemp.fnd TSTART=0 TSTOP=100 NEVAL=200resulting in the FINEDATA control stream file nmtemp.fnd:
$INPUT C SET ID JID TIME DV=CONC AMT=DOSE RATE EVID MDV CMT CLX V1X QX V2X
SDIX SDSX
$DATA nmtemp.csv IGNORE=C
$FINEDATA AXIS=TIME(LIN) TSTOP=100 TSTART=0 NEVAL=200 FILE=nmtemp2.csvNote that only words that match the variable list at the nmtemplate command line, and have enclosing brackets <>, will be replaced with the suggested values. The values may also be text with no spaces in them.
These two scripts could be combined to provide a means of creating individual simulated curves. Consider the following DOS patch script (which could also be converted to an R/S-PLUS script or function), nmtemp.bat:
nmtemplate.exe nmtemp.fnt nmtemp.fnd TSTART=%1 TSTOP=%2 NEVAL=%3
finedata.exe nmtemp.fnd
nmtemplate.exe nmtemp.nmt nmtemp.ctl NMID=%4
$nmfe75.bat nmtemp.ctl nmtemp.res -prdefaultWhere %1 through %4 are the DOS command line substitution parameters. So the script could be executed as follows:
Call nmtemp.bat 0 100 200 34Then, a program such as R, S-PLUS, or S-ADAPT, can read in the results from nmtemp.tab, and plot them.
Another feature of nmtemplate is that the user may request a random number to be generated to serve as a value, by referring to ~R(a1,a2,a3). R(a1,a2,a3) is a special function of nmtemplate, which obtains a uniform random variate between a1 and a2. If a seed a3 is given that is not 0, it means to initialize the seed. The initialization should be done once in a series. For example:
The following line sets the seed:
nmtemplate wexample12.nmt dummy.ctl SAMPLE=~R(1,10000,113345)with a throw-away result file dummy.ctl. Then one could perform a for loop in a DOS batch file to generate a series of control stream files with different starting seeds:
for /l %%n in (1,1,9) do nmtemplate wexample12.nmt wexample12_%%n.ctl SAMPLE=~R(%%n000,%%n999,0)where for /l %%n in (1,1,9) is a DOS command generating n starting at 1, incrementing by 1, and ending at 9. When n=3, for example, ~R(%%n000,%%n999,0) will be ~R(3000,3999,0), generating a random number between 3000 and 3999, to be substituted wherever <SAMPLE> shows up in the template file wexample12.nmt.
The template file wexample12.nmt may contain:
$EST METHOD=CHAIN FILE=wexample12.txt NSAMPLE=0 ISAMPLE=<SAMPLE>and the resulting files wexample12_1.ctl through wexample12_9.ctl will contain random ISAMPLE values, such as:
wexample12_1.ctl:
$EST METHOD=CHAIN FILE=wexample12.txt NSAMPLE=0 ISAMPLE=1345wexample12_2.ctl:
$EST METHOD=CHAIN FILE=wexample12.txt NSAMPLE=0 ISAMPLE=2456wexample12_3.ctl:
$EST METHOD=CHAIN FILE=wexample12.txt NSAMPLE=0 ISAMPLE=3089etc. It should be pointed out that this example, in which nmtemplate is used to create a random variable for substitution into ISAMPLE, can easily be done in NM73 using the ISAMPEND and SELECT=3 options for $EST METHOD=CHAIN or $CHAIN (see I.66 Method for creating several instances for a problem starting at different randomized initial positions: $EST METHOD=CHAIN and $CHAIN Records).
You may use STDIN for console input, and STDOUT for console output, useful for piping commands together, for example:
nmtemplate doetest.ctl STDOUT IMAX=3 JMAX=4 KMAX=14 | doexpand >doetestq.ctl22 NEFF and NEFFI Utility Programs (NM74)
In MCMC sampling, often, the number of effective samples is desired to be determined, which takes into account the correlation between samples. You may collect 10000 samples during the stationary phase, but these may be only equivalent to 100 independent samples. The NEFF utility analyzes EXT files containing Bayesian samples generated from NONMEM. Samples from all tables in one or more files can be collected, and overall NEFF statistics will be generated. For each file, enter the minimum and maximum iteration number value. When finished entering file names, enter the word STOP (in uppercase). For example:
c:\nm74a6\util>neff
ENTER EXT FILE name, MIN, MAX
stanrb_169.ext,0,10000
stanrb_169.ext 0 10000
10000
ENTER EXT FILE name, MIN, MAX
stanrb_171.ext,0,10000
stanrb_171.ext 0
ENTER EXT FILE name, MIN, MAX
STOP
PARAMETER MEAN STD N NEFF1 %NEFF1 NEFF2 %NEFF2 RHAT
THETA1 -2.505708E-03 4.006207E-03 4002 4002.0 100.00 4002.0 100.00 1.000
THETA2 2.921375E-02 4002 2655.4 66.35 2801.9 70.01 1.000
THETA3 3.68305 4.326550E-02 4002 368.8 9.22 406.7 10.16 1.006
THETA4 -5.01367 0.142392 4002 1843.7 46.07 1866.7 46.65 1.000
THETA5 -0.988017 0.162227 4002 1843.5 46.06 1885.3 47.11 1.000
SIGMA(1,1) -1.12213 0.265440 4002 4002.0 100.00 4002.0 100.00 1.000
OMEGA(1,1) 6.782039E-04 4002 2274.6 56.84 2492.6 62.28 1.000
OMEGA(2,1) 15.8272 2.568944E-03 4002 303.5 7.58 512.5 12.81 1.009
OMEGA(2,2) 1.028170E-02 1.291445E-02 4002 296.2 7.40 558.0 13.94 1.010
OMEGA(3,1) 2.716656E-02 1.690585E-03 4002 226.1 5.65 296.1 7.40 1.007
OMEGA(3,2) 0.175062 5.949785E-03 4002 7.7 0.19 153.9 3.85 1.087
OMEGA(3,3) 7.886005E-04 5.804516E-03 4002 5.2 0.13 149.4 3.73 1.134
MCMCOBJ -4.639637E-02 4002 2.1 0.05 76.6 1.91 1.481
4.143005E-02 473.698
26650.1You may add a delimiter to the command line of neff, if for example, the input files are comma delimited:
neff ,The NEFF1 statistics are based on the method of evaluation in the appendices of the STAN reference manual [20]. The NEFF2 statistics are based on Appendix A of [19].
22.1 NEFFI for individual parameters
If you wish to analyze phi or eta values, you can select the BAYES_PHI_STORE=1 option for $EST (see BAYES_PHI_STORE=0 (default) (NM75) in section I.41 Full Markov Chain Monte Carlo (MCMC) Bayesian Analysis Method), which will produce a root.iph file that the NEFFI utility can read.
For outputting other individual parameters, you will need to add additional lines in $PK or $PRED, using the BAYES_EXTRA_REQUEST signal (see section I.70 Obtaining random samples of individual predicted values and individual parameters during MCMC Bayesian Analysis, and example8 as well, near the end of this document), and output the desired individual parameters. For example:
$PRED
include nonmem_reserved_general
BAYES_EXTRA_REQUEST=1
…
(code here)
…
PHI1=ETA(1)+MU_1
PHI2=ETA(2)+MU_2
PHI3=ETA(3)+MU_3
IF(ICALL==1) THEN
" OPEN(UNIT=50,FILE='stanb_new296.eta')
" WRITE(50,'(A)') 'TABLE NO. 2: NUTS Bayesian Analysis: Phi Values'
" WRITE(50,'(A12,1X,A14,3(1X,A12))') 'ITERATION','ID','PHI1','PHI2','PHI3'
ENDIF
; GET THE FIRST RECORD OF EACH SUBJECT, AND ONLY WHEN IBMETHOD=EST_BAYES
; AND BAYES_METHOD=1 (NUTS),
; AND DURING THE STATIONARY ITERATIONS (ITER_REPORT>0)
IF(BAYES_EXTRA==1 .AND.IBMETHOD==EST_BAYES.AND.BAYES_METHOD==1. &
AND.NEWIND/=2.AND.ITER_REPORT>=0) THEN
" WRITE(50,'(I12,1X,F14.0,3(1X,1PG12.5))') ITER_REPORT,ID,PHI1,PHI2,PHI3
ENDIFThis will produce a file of the following structure:
TABLE NO. 2: NUTS Bayesian Analysis: Phi Values
ITERATION ID PHI1 PHI2 PHI3
-4.4562
0 1. 4.55347E-02 -0.10868 -4.4000
-4.2769
0 2. 3.51615E-02 -0.38840 -4.3652
-4.3818
0 3. 8.05289E-02 -0.18820 -4.2215
-4.3358
0 4. -0.10113 -1.1717 -4.2898
-4.5667
0 5. -6.03206E-02 -0.58953 -4.3962
-4.1605
0 6. -9.53758E-02 -1.3678 -4.5920
-4.2070
0 7. -6.60832E-02 -0.67833
0 8. -7.32508E-02 -1.0046
0 9. 2.12179E-02 -0.41888
0 10. 0.13819 -0.20922
0 11. 0.16236 -0.81138
0 12. 5.78395E-02 -0.30828
0 13. 5.00854E-02 -1.0830Notice that there needs to be a table header with "TABLE NO. " specified, followed by the column header, the first two items begin iteration, then ID (or subject number NIREC), followed by a list of desired items, in accordance with the write statement in the control stream file.
A file of this structure can be read by the utility program NEFFI (for NEFF individual), for example (a list of table files may be given as with NEFF, ending with STOP):
neffi
ENTER EXT FILE name, MIN, MAX
stanb_new296.eta,0,10000
STOPThe efficiency and correlation statistics will be printed out, with labels of type XXXX_itemname, where XXXX is the subject number.