Design evaluation and optimization in discrete space

Overview

This example is based on a landmark PK/PD population study of Tolmetin (a non-steroidal anti-inflammatory drug) in rats (Flores-Murrieta et al. 1998).

The model consists of:

Original experimental design

The original study involved 6 parallel groups of rats (n ≥ 6 per group) receiving single oral doses of 1, 3.2, 10, 31.6, 56.2, or 100 mg/kg per os. Blood sampling and drug response (DI score) evaluation were conducted at 0, 15, 30, and 45 min and at 1, 1.25, 1.5, 2, 3 and 4 hours after administration (nine non-zero sampling times in hours: 0.25–4 h).

Objectives

  1. Evaluation: compute the Population Fisher Information Matrix (FIM) for the original design to assess the precision of structural parameters and variance components (RSE%, Shrinkage). A Bayesian FIM is also computed for comparison.
  2. Optimization: identify the D-optimal design for a total of 30 rats by selecting dose levels and sampling times from a discrete candidate set. The Fedorov-Wynn and Multiplicative algorithms are compared.

Optimisation results are computed by example01_execute.R (run once, then cached as .RDS in data/). HTML reports are written to results/. Reports are also available at https://github.com/packagePFIM

Design evaluation

Model equations

The PKPD model is defined as a system of Ordinary Differential Equations (ODEs) using named character strings. PFIM performs symbolic differentiation on these strings to derive the sensitivity equations required for FIM computation.

Convention:

Equation 1 (PK) — one-compartment model with first-order oral absorption:

\[\frac{dC_c}{dt} = \frac{\mathrm{dose_{RespPK}}}{V} \cdot k_a \cdot e^{-k_a t} - \frac{Cl}{V} \cdot C_c\]

Symbol Description
V volume of distribution (L)
ka first-order absorption rate constant (h⁻¹)
Cl total clearance (L/h)
Cc plasma drug concentration — state variable (mcg/mL)

Equation 2 (PD) — indirect response model, inhibition of production (Type I):

\[\frac{dE}{dt} = R_{in} \left(1 - I_{max} \frac{C_c^\gamma}{C_c^\gamma + IC_{50}^\gamma}\right) - k_{out} \cdot E\]

Symbol Description
Rin baseline production rate of the inflammation score (h⁻¹)
Imax maximum fractional inhibition (dimensionless, 0 < Imax ≤ 1)
IC50 concentration producing 50% of Imax (mcg/mL)
gamma Hill coefficient (dimensionless)
kout first-order elimination rate of the effect (h⁻¹)
E pharmacodynamic response (DI inflammation score) — state variable

Steady-state note: at t = 0 (Cc = 0) the system is at equilibrium: E(0) = Rin / kout (here 614 / 6.14 = 100).


modelEquations = list(
  "Deriv_Cc" = "dose_RespPK/V*ka*exp(-ka*t) - Cl/V*Cc",
  "Deriv_E"  = "Rin*(1-Imax*(Cc**gamma)/(Cc**gamma + IC50**gamma))-kout*E"
)

Model parameters

Parameters are specified via their population typical value (fixed effect, mu) and inter-individual variability (IIV, omega). PFIM assumes a log-normal distribution for all parameters, guaranteeing strict positivity.

The parameter vector estimated by the population FIM is:

\[\theta = \{\mu_V, \mu_{Cl}, \mu_{kout}, \mu_{Imax}, \mu_{IC50}, \mu_{gamma}, \omega^2_V, \omega^2_{Cl}, \omega^2_{kout}, \omega^2_{Imax}, \omega^2_{IC50}, \omega^2_{gamma}\}\]

Parameters with fixedMu = TRUE are considered known constants and are excluded from the FIM. Parameters with omega = 0 have no IIV component and their variance is not estimated.

Name Description mu omega fixedMu
V Volume of distribution (L) 0.74 0.316 FALSE
Cl Total clearance (L/h) 0.28 0.456 FALSE
ka Absorption rate constant (h⁻¹) 10 0 TRUE
kout Effect elimination rate (h⁻¹) 6.14 0.947 FALSE
Rin Baseline production rate (h⁻¹) 614 0 TRUE
Imax Maximum inhibition (-) 0.76 0.439 FALSE
IC50 Potency (mcg/mL) 9.22 0.452 FALSE
gamma Hill coefficient (-) 2.77 1.761 FALSE

Rationale for fixed parameters:


modelParameters = list(

  ModelParameter(name         = "V",
                 distribution = LogNormal(mu = 0.74,  omega = 0.316)),

  ModelParameter(name         = "Cl",
                 distribution = LogNormal(mu = 0.28,  omega = 0.456)),

  ModelParameter(name         = "ka",
                 distribution = LogNormal(mu = 10,    omega = 0),
                 fixedMu      = TRUE),

  ModelParameter(name         = "kout",
                 distribution = LogNormal(mu = 6.14,  omega = 0.947)),

  ModelParameter(name         = "Rin",
                 distribution = LogNormal(mu = 614,   omega = 0),
                 fixedMu      = TRUE),

  ModelParameter(name         = "Imax",
                 distribution = LogNormal(mu = 0.76,  omega = 0.439)),

  ModelParameter(name         = "IC50",
                 distribution = LogNormal(mu = 9.22,  omega = 0.452)),

  ModelParameter(name         = "gamma",
                 distribution = LogNormal(mu = 2.77,  omega = 1.761))
)

Residual error models

Two distinct residual error models are specified, one per response.

Combined1(output, sigmaInter, sigmaSlope) — combined additive + proportional model:

\[\mathrm{SD}(\epsilon) = \sigma_{inter} + \sigma_{slope} \cdot f(\theta, \xi)\]

Setting sigmaInter = 0 reduces it to a pure proportional model: SD(ε_PK) = 0.21 × Cc (21% proportional error). This is appropriate for plasma concentrations where measurement error scales with the signal magnitude across the dynamic range.

Constant(output, sigmaInter) — additive error model: SD(ε_PD) = 9.6 DI units. This is appropriate for bounded inflammation scores whose measurement precision does not depend on the response level.


errorModelRespPK = Combined1(output = "RespPK", sigmaInter = 0,   sigmaSlope = 0.21)
errorModelRespPD = Constant( output = "RespPD", sigmaInter = 9.6)

modelError = list(errorModelRespPK, errorModelRespPD)

Sampling times

Nine observation times (hours) are used for both responses, spanning:

Time 0 is omitted: Cc(0) = 0 and E(0) = 100 are fixed initial conditions that carry no information about model parameters (zero sensitivity).


samplingTimesRespPK = SamplingTimes(
  outcome   = "RespPK",
  samplings = c(0.25, 0.5, 0.75, 1, 1.25, 1.5, 2, 3, 4)
)

samplingTimesRespPD = SamplingTimes(
  outcome   = "RespPD",
  samplings = c(0.25, 0.5, 0.75, 1, 1.25, 1.5, 2, 3, 4)
)

Arms (dose groups)

Six arms correspond to the original dose levels, converted from mg/kg to absolute doses for a 200 g rat (dose_mg = dose_mg/kg × 0.200 kg).

Arm mg/kg Absolute dose Subjects
1 1.0 0.20 mg 6
2 3.2 0.64 mg 6
3 10.0 2.00 mg 6
4 31.6 6.32 mg 6
5 56.2 11.24 mg 6
6 100.0 20.00 mg 6

Design summary: 6 arms × 6 subjects = 36 subjects total.

Initial conditions:


administrationRespPK1 = Administration(
  outcome  = "RespPK",
  timeDose = c(0),
  dose     = c(0.2)
)

arm1 = Arm(
  name             = "0.2mg Arm",
  size             = 6,
  administrations  = list(administrationRespPK1),
  samplingTimes    = list(samplingTimesRespPK, samplingTimesRespPD),
  initialCondition = list("Cc" = 0, "E" = 100)
)

administrationRespPK2 = Administration(outcome = "RespPK", timeDose = c(0), dose = c(0.64))

arm2 = Arm(
  name             = "0.64mg Arm",
  size             = 6,
  administrations  = list(administrationRespPK2),
  samplingTimes    = list(samplingTimesRespPK, samplingTimesRespPD),
  initialCondition = list("Cc" = 0, "E" = 100)
)

administrationRespPK3 = Administration(outcome = "RespPK", timeDose = c(0), dose = c(2))

arm3 = Arm(
  name             = "2mg Arm",
  size             = 6,
  administrations  = list(administrationRespPK3),
  samplingTimes    = list(samplingTimesRespPK, samplingTimesRespPD),
  initialCondition = list("Cc" = 0, "E" = 100)
)

administrationRespPK4 = Administration(outcome = "RespPK", timeDose = c(0), dose = c(6.32))

arm4 = Arm(
  name             = "6.32mg Arm",
  size             = 6,
  administrations  = list(administrationRespPK4),
  samplingTimes    = list(samplingTimesRespPK, samplingTimesRespPD),
  initialCondition = list("Cc" = 0, "E" = 100)
)

administrationRespPK5 = Administration(outcome = "RespPK", timeDose = c(0), dose = c(11.24))

arm5 = Arm(
  name             = "11.24mg Arm",
  size             = 6,
  administrations  = list(administrationRespPK5),
  samplingTimes    = list(samplingTimesRespPK, samplingTimesRespPD),
  initialCondition = list("Cc" = 0, "E" = 100)
)

administrationRespPK6 = Administration(outcome = "RespPK", timeDose = c(0), dose = c(20))

arm6 = Arm(
  name             = "20mg Arm",
  size             = 6,
  administrations  = list(administrationRespPK6),
  samplingTimes    = list(samplingTimesRespPK, samplingTimesRespPD),
  initialCondition = list("Cc" = 0, "E" = 100)
)

Design assembly

Design() aggregates all arms into a single experimental design object.


design1 = Design(
  name = "design1",
  arms = list(arm1, arm2, arm3, arm4, arm5, arm6)
)

Population and Bayesian FIM evaluation

The Evaluation() constructor specifies the full statistical model:

Argument Role
modelEquations user-defined ODE system
modelParameters fixed effects + IIV
modelError intra-individual error
outputs named list mapping outcome labels to ODE state variables ("RespPK" → Cc, "RespPD" → E)
designs list of Design objects to evaluate
fimType "population" estimates θ = {mu, ω²}; "Bayesian" is the individual FIM regularized by a prior on ω²
odeSolverParameters passed to deSolve::lsoda; tight tolerances (1e-8) are required because the PD sub-model is moderately stiff (kout = 6.14 h⁻¹ implies rapid equilibration of the effect)

PFIM applies a First-Order (FO) linearization, expanding the individual model in a first-order Taylor series around the typical values mu. The population FIM has dimension p × p, where p = number of estimable parameters (here p = 12: 6 fixed effects + 6 variance components).

Population FIM


evaluationPop = Evaluation(
  name                = "evaluation",
  modelEquations      = modelEquations,
  modelParameters     = modelParameters,
  modelError          = modelError,
  outputs             = list("RespPK" = "Cc", "RespPD" = "E"),
  designs             = list(design1),
  fimType             = "population",
  odeSolverParameters = list(atol = 1e-8, rtol = 1e-8)
)

evaluationPop = run(evaluationPop)

Bayesian FIM

The Bayesian FIM augments the individual FIM with the inverse prior covariance matrix (Ω⁻¹), acting as a regularization term. This is equivalent to a Maximum A Posteriori (MAP) estimation framework and is relevant when prior information on ω² is available. All other arguments are identical to the population evaluation above.


evaluationBay = Evaluation(
  name                = "evaluation",
  modelEquations      = modelEquations,
  modelParameters     = modelParameters,
  modelError          = modelError,
  outputs             = list("RespPK" = "Cc", "RespPD" = "E"),
  designs             = list(design1),
  fimType             = "Bayesian",
  odeSolverParameters = list(atol = 1e-8, rtol = 1e-8)
)

evaluationBay = run(evaluationBay)

Results: display and export

Accessor functions for retrieving results from Evaluation / Optimization objects:

Function Description
show(x) formatted summary of all statistical metrics
getFisherMatrix(x) retrieves the FIM
getCorrelationMatrix(x) normalized FIM to identify parameter collinearity
getSE(x) asymptotic Standard Errors (SE)
getRSE(x) Relative Standard Errors (RSE%)
getShrinkage(x) Shrinkage (%) for random effects
getDeterminant(x) determinant of the FIM
getDcriterion(x) D-optimality criterion of the FIM 

show(evaluationPop)
writeLines(capture.output(show(evaluationPop)),
           file.path(paths$outputs, "vignette1_evaluation_populationFIM_show.txt"))
fisherMatrix = getFisherMatrix(evaluationPop)
getCorrelationMatrix(evaluationPop)
getSE(evaluationPop)
getRSE(evaluationPop)
getShrinkage(evaluationPop)
getDeterminant(evaluationPop)
getDcriterion(evaluationPop)

show(evaluationBay)
writeLines(capture.output(show(evaluationBay)),
           file.path(paths$outputs, "vignette1_evaluation_BayesianFIM_show.txt"))
fisherMatrix = getFisherMatrix(evaluationBay)
getCorrelationMatrix(evaluationBay)
getSE(evaluationBay)
getRSE(evaluationBay)
getShrinkage(evaluationBay)
getDeterminant(evaluationBay)
getDcriterion(evaluationBay)
*************************************** 
  Population Fisher Matrix 
*************************************** 

                        μ_V          μ_Cl      μ_kout      μ_Imax      μ_IC50     μ_gamma         ω²_V        ω²_Cl      ω²_kout      ω²_Imax    ω²_IC50   ω²_gamma σ_slope_RespPK σ_inter_RespPD
μ_V            583.89694102    9.19036099 -0.02315464 -3.35612090  0.32079096 -0.24507445 0.000000e+00 0.000000e+00 0.000000e+00  0.000000000  0.0000000 0.00000000   0.000000e+00    0.000000000
μ_Cl             9.19036099 2110.23408215 -0.01715815 -0.03198433  0.52533166  0.35319017 0.000000e+00 0.000000e+00 0.000000e+00  0.000000000  0.0000000 0.00000000   0.000000e+00    0.000000000
μ_kout          -0.02315464   -0.01715815  1.04306085  0.63336569 -0.03141002 -0.04209541 0.000000e+00 0.000000e+00 0.000000e+00  0.000000000  0.0000000 0.00000000   0.000000e+00    0.000000000
μ_Imax          -3.35612090   -0.03198433  0.63336569 60.63827631 -2.70317529  1.54206818 0.000000e+00 0.000000e+00 0.000000e+00  0.000000000  0.0000000 0.00000000   0.000000e+00    0.000000000
μ_IC50           0.32079096    0.52533166 -0.03141002 -2.70317529  0.65184243  0.07044445 0.000000e+00 0.000000e+00 0.000000e+00  0.000000000  0.0000000 0.00000000   0.000000e+00    0.000000000
μ_gamma         -0.24507445    0.35319017 -0.04209541  1.54206818  0.07044445  0.74226627 0.000000e+00 0.000000e+00 0.000000e+00  0.000000000  0.0000000 0.00000000   0.000000e+00    0.000000000
ω²_V             0.00000000    0.00000000  0.00000000  0.00000000  0.00000000  0.00000000 1.419936e+03 5.071130e-02 5.312959e-04  0.174193007  0.4767726 0.01121281   1.694963e+02    0.022484025
ω²_Cl            0.00000000    0.00000000  0.00000000  0.00000000  0.00000000  0.00000000 5.071130e-02 3.801555e+02 9.340249e-05  0.003376237  0.1153787 0.00215635   3.442924e+01    0.004590882
ω²_kout          0.00000000    0.00000000  0.00000000  0.00000000  0.00000000  0.00000000 5.312959e-04 9.340249e-05 2.149292e+01  0.386744465  0.1693550 0.01150713   8.462700e-03    0.050138553
ω²_Imax          0.00000000    0.00000000  0.00000000  0.00000000  0.00000000  0.00000000 1.741930e-01 3.376237e-03 3.867445e-01 41.669312017 16.1932689 0.55767351   1.413532e+00    0.818888973
ω²_IC50          0.00000000    0.00000000  0.00000000  0.00000000  0.00000000  0.00000000 4.767726e-01 1.153787e-01 1.693550e-01 16.193268935 83.8259799 0.15388072   8.155428e+00    1.229640291
ω²_gamma         0.00000000    0.00000000  0.00000000  0.00000000  0.00000000  0.00000000 1.121281e-02 2.156350e-03 1.150713e-02  0.557673507  0.1538807 0.69850714   1.549203e-01    0.107175844
σ_slope_RespPK   0.00000000    0.00000000  0.00000000  0.00000000  0.00000000  0.00000000 1.694963e+02 3.442924e+01 8.462700e-03  1.413531724  8.1554276 0.15492032   1.146177e+04    0.328767259
σ_inter_RespPD   0.00000000    0.00000000  0.00000000  0.00000000  0.00000000  0.00000000 2.248402e-02 4.590882e-03 5.013855e-02  0.818888973  1.2296403 0.10717584   3.287673e-01    5.297969806

*************************************** 
  Fixed effects (μ) 
*************************************** 

                 μ_V          μ_Cl      μ_kout      μ_Imax      μ_IC50     μ_gamma
μ_V     583.89694102    9.19036099 -0.02315464 -3.35612090  0.32079096 -0.24507445
μ_Cl      9.19036099 2110.23408215 -0.01715815 -0.03198433  0.52533166  0.35319017
μ_kout   -0.02315464   -0.01715815  1.04306085  0.63336569 -0.03141002 -0.04209541
μ_Imax   -3.35612090   -0.03198433  0.63336569 60.63827631 -2.70317529  1.54206818
μ_IC50    0.32079096    0.52533166 -0.03141002 -2.70317529  0.65184243  0.07044445
μ_gamma  -0.24507445    0.35319017 -0.04209541  1.54206818  0.07044445  0.74226627

*************************************** 
  Variance components (ω², γ², σ) 
*************************************** 

                       ω²_V        ω²_Cl      ω²_kout      ω²_Imax    ω²_IC50   ω²_gamma σ_slope_RespPK σ_inter_RespPD
ω²_V           1.419936e+03 5.071130e-02 5.312959e-04  0.174193007  0.4767726 0.01121281   1.694963e+02    0.022484025
ω²_Cl          5.071130e-02 3.801555e+02 9.340249e-05  0.003376237  0.1153787 0.00215635   3.442924e+01    0.004590882
ω²_kout        5.312959e-04 9.340249e-05 2.149292e+01  0.386744465  0.1693550 0.01150713   8.462700e-03    0.050138553
ω²_Imax        1.741930e-01 3.376237e-03 3.867445e-01 41.669312017 16.1932689 0.55767351   1.413532e+00    0.818888973
ω²_IC50        4.767726e-01 1.153787e-01 1.693550e-01 16.193268935 83.8259799 0.15388072   8.155428e+00    1.229640291
ω²_gamma       1.121281e-02 2.156350e-03 1.150713e-02  0.557673507  0.1538807 0.69850714   1.549203e-01    0.107175844
σ_slope_RespPK 1.694963e+02 3.442924e+01 8.462700e-03  1.413531724  8.1554276 0.15492032   1.146177e+04    0.328767259
σ_inter_RespPD 2.248402e-02 4.590882e-03 5.013855e-02  0.818888973  1.2296403 0.10717584   3.287673e-01    5.297969806

********************************************* 
  Determinant, condition numbers and D-criterion 
 *********************************************** 

Determinant: 4.246897e+22 
D-criterion: 41.33242 
Condition number (fixed effects): 4677.284 
Condition number (variance components): 16644.23 

*************************************** 
  Parameters estimation 
*************************************** 

Parameter               Value           SE     RSE(%)
μ_V                 0.740000  0.041396101   5.594068
μ_Cl                0.280000  0.021772569   7.775917
μ_kout              6.140000  0.984617978  16.036123
μ_Imax              0.760000  0.149908136  19.724755
μ_IC50              9.220000  1.409222187  15.284406
μ_gamma             2.770000  1.227830929  44.326026
ω²_V               0.099856  0.026561314  26.599617
ω²_Cl              0.207936  0.051295421  24.668851
ω²_kout            0.896809  0.215721035  24.054290
ω²_Imax            0.192721  0.162032545  84.076227
ω²_IC50            0.204304  0.113693973  55.649411
ω²_gamma           3.101121  1.204551613  38.842458
σ_slope_RespPK       0.210000  0.009350449   4.452595
σ_inter_RespPD       9.600000  0.436125367   4.542973

[1] 41.33242
*************************************** 
 Bayesian Fisher Matrix 
*************************************** 

               μ_V        μ_Cl    μ_kout      μ_Imax     μ_IC50     μ_gamma
μ_V     89.9225529  10.8651606 -2.304353 -0.99425042  1.1844166 -0.14172892
μ_Cl    10.8651606 110.2098299 -3.726483 -0.48861216  1.7730592  0.46195110
μ_kout  -2.3043531  -3.7264832 69.674268  6.31774142 -6.0308411 -2.35011313
μ_Imax  -0.9942504  -0.4886122  6.317741  7.15950553 -1.4828623 -0.00546153
μ_IC50   1.1844166   1.7730592 -6.030841 -1.48286233  7.8521422  0.32022113
μ_gamma -0.1417289   0.4619511 -2.350113 -0.00546153  0.3202211  0.72337746

*************************************** 
 Fixed effects 
*************************************** 

               μ_V        μ_Cl    μ_kout      μ_Imax     μ_IC50     μ_gamma
μ_V     89.9225529  10.8651606 -2.304353 -0.99425042  1.1844166 -0.14172892
μ_Cl    10.8651606 110.2098299 -3.726483 -0.48861216  1.7730592  0.46195110
μ_kout  -2.3043531  -3.7264832 69.674268  6.31774142 -6.0308411 -2.35011313
μ_Imax  -0.9942504  -0.4886122  6.317741  7.15950553 -1.4828623 -0.00546153
μ_IC50   1.1844166   1.7730592 -6.030841 -1.48286233  7.8521422  0.32022113
μ_gamma -0.1417289   0.4619511 -2.350113 -0.00546153  0.3202211  0.72337746

*********************************************** 
 Determinant, condition numbers and D-criterion 
*********************************************** 

Determinant: 20333648.444364 
D-criterion: 16.5209830299747 
Condition number of the fixed effects: 182.897110733334 

*************************************** 
 Shrinkage 
*************************************** 

               μ_V     μ_Cl   μ_kout   μ_Imax   μ_IC50 μ_gamma
Shrinkage 11.31068 4.440497 2.041112 81.25003 68.55303 50.9391

*************************************** 
 Parameters estimation 
*************************************** 

Parameter             Value           SE     RSE(%)
μ_V                   0.74   0.07864357  10.627510
μ_Cl                  0.28   0.02690535   9.609054
μ_kout                6.14   0.83071446  13.529551
μ_Imax                0.76   0.30073909  39.570933
μ_IC50                9.22   3.45050512  37.424134
μ_gamma               2.77   3.48148698 125.685451

[1] 16.52098

Diagnostic plots

plotEvaluation(), plotSensitivityIndices(), plotSE(), and plotRSE() are the PFIM entry points for evaluation graphics:

plotOptions controls axis labels in all PFIM graphics.


plotOptions = list(unitTime = c("hour"), unitOutcomes = c("mcg/mL", "DI%"))

plotsEval1_eval = plotEvaluation(evaluationPop, plotOptions)
plotsEval1_si   = plotSensitivityIndices(evaluationPop, plotOptions)

plotOutcomesEvaluationRespPK = plotsEval1_eval$design1$`20mg Arm`$RespPK
plotOutcomesEvaluationRespPD = plotsEval1_eval$design1$`20mg Arm`$RespPD
plotSensitivityIndice_RespPK_Cl = plotsEval1_si$design1$`20mg Arm`$RespPK$Cl
plotSensitivityIndice_RespPK_V = plotsEval1_si$design1$`20mg Arm`$RespPK$V
plotEval_SE = PFIM::plotSE(evaluationPop)
plotEval_RSE = PFIM::plotRSE(evaluationPop)

ggsave(file.path(paths$figures, "vignette1_evaluation_populationFim_design1_arm20mg_RespPK.pdf"),
       plotOutcomesEvaluationRespPK, width = 8, height = 5)

HTML report


outputFile = "Example01_EvaluationPopFIM.html"
Report(evaluationPop, paths$reports, outputFile, plotOptions)

Design optimization

Objectives and constraints

Building on the evaluation above, we now seek an optimal design for a future study under practical constraints:

Algorithms compared

Both algorithms operate in the discrete candidate space and maximize the D-criterion of the population FIM.

Fedorov-Wynn (FW): an exact exchange algorithm that iteratively adds the elementary protocol (dose × sampling schedule pair) that most increases the FIM determinant, then removes the least contributing one. It converges to a D-optimal design on the discrete support of candidate protocols.

Multiplicative Algorithm (MA): a continuous relaxation method that assigns and iteratively updates weights to all candidate protocols. Weights below a threshold are zeroed at convergence, yielding a sparse approximate D-optimal design. Unlike FW, MA does not require initial elementary protocols and explores the full candidate space simultaneously.

Runtime note

Both algorithms require approximately 10–15 minutes to run on a standard workstation. During vignette rendering, example01_execute.R runs each optimization once (with showProcess = FALSE) and saves the result to data/; subsequent renders load the cached .RDS files.

Initial administration and candidate sampling grids

Starting dose for the constrained arm. The optimizer will reassign doses to subjects from the discrete set defined in the administration constraints below.


administrationRespPK = Administration(outcome = "RespPK", timeDose = c(0), dose = c(6.32))

These are the full sets of candidate time points (hours) from which the optimizer will select the most informative subset, subject to the constraints defined in the next sections.


samplingTimesRespPK = SamplingTimes(
  outcome   = "RespPK",
  samplings = c(0.25, 0.75, 1, 1.5, 2, 4, 6)
)

samplingTimesRespPD = SamplingTimes(
  outcome   = "RespPD",
  samplings = c(0.25, 0.75, 1.5, 2, 3, 6, 8, 12)
)

Sampling time constraints

SamplingTimeConstraints(outcome, initialSamplings, fixedTimes, numberOfsamplingsOptimisable, ...) restricts which time points can be assigned; one object per outcome in the arm.

RespPK constraints:

RespPD constraints:


samplingConstraintsRespPK = SamplingTimeConstraints(
  outcome                      = "RespPK",
  initialSamplings             = c(0.25, 0.75, 1, 1.5, 2, 4, 6),
  fixedTimes                   = c(0.25, 4),
  numberOfsamplingsOptimisable = 4
)

samplingConstraintsRespPD = SamplingTimeConstraints(
  outcome                      = "RespPD",
  initialSamplings             = c(0.25, 0.75, 1.5, 2, 3, 6, 8, 12),
  fixedTimes                   = c(2, 6),
  numberOfsamplingsOptimisable = 4
)

Initial elementary protocols for Fedorov-Wynn

An elementary protocol is a (dose, sampling schedule) combination for a homogeneous subgroup of subjects. In a PK/PD arm the sampling schedule is multi-outcome: PK times and PD times are concatenated on the FW candidate grid.

The initial support used here is one protocol (with proportionsOfSubjects of length 1):

Pass it as a nested list (one element = one support point), or equivalently as a single flat vector of length 8. A vignette-style list(pk, pd) with one proportion is also accepted.

initialElementaryProtocols is passed to optimizerParameters$elementaryProtocols inside the Optimization() call; it is only used by FedorovWynnAlgorithm.


initialElementaryProtocols = list(
  list(
    c(0.25, 0.75, 1, 4),
    c(1.5, 2, 6, 12)
  )
)

Dose constraints

AdministrationConstraints(outcome, doses) restricts which dose values the optimizer can assign to each elementary protocol. The discrete set (in mg) corresponds to the 6 dose levels of the original study.


administrationConstraintsRespPK = AdministrationConstraints(
  outcome = "RespPK",
  doses   = list(0.2, 0.64, 2, 6.32, 11.24, 20)
)

Constrained arm and design

The arm encodes all constraints simultaneously and serves as the template that both optimization algorithms will operate on.

The initial condition for E is specified as "Rin/kout" (a formula string) rather than the numeric value 100. PFIM evaluates this expression at the typical parameter values, giving Eâ‚€ = 614 / 6.14 = 100. This ensures that the initial condition remains consistent with the model structure if parameter estimates are updated.

numberOfArms is the upper bound on the number of distinct elementary protocols the optimizer can create.


armConstraint = Arm(
  name                       = "armConstraint",
  size                       = 30,
  administrations            = list(administrationRespPK),
  samplingTimes              = list(samplingTimesRespPK, samplingTimesRespPD),
  administrationsConstraints = list(administrationConstraintsRespPK),
  samplingTimesConstraints   = list(samplingConstraintsRespPK, samplingConstraintsRespPD),
  initialCondition           = list("Cc" = 0, "E" = "Rin/kout")
)

designConstraint = Design(
  name         = "designConstraint",
  arms         = list(armConstraint),
  numberOfArms = 30
)

numberOfSubjects      = c(30)
proportionsOfSubjects = c(30) / 30

For very large dose × sampling grids, cap FIM evaluations (deterministic subsample):


pfim_set_option(constraints.maxTasks = 500)

Fedorov-Wynn algorithm

The Fedorov-Wynn algorithm is an exact combinatorial exchange method for finding D-optimal designs in a discrete candidate space. It proceeds as:

  1. Start from the user-supplied initial elementary protocols.
  2. Evaluate the directional derivative of the D-criterion for every candidate protocol not currently in the design support.
  3. Add the protocol with the highest derivative (greedy ascent step).
  4. Remove the protocol that contributes least to the current FIM.
  5. Repeat steps 2–4 until the improvement in |FIM| falls below a tolerance.

FedorovWynnAlgorithm optimizerParameters:

Parameter Description
elementaryProtocols list of protocols; each protocol is a flat numeric vector (full grid row) or a list of per-outcome vectors (e.g. PK then PD)
numberOfSubjects integer vector; total N to distribute across protocols
proportionsOfSubjects numeric vector summing to 1; initial allocation fractions
showProcess logical; if TRUE, prints per-iteration D-criterion values 

optimizationFWPopFIM = Optimization(
  name                = "PKPD_ODE_multi_doses_populationFIM",
  modelEquations      = modelEquations,
  modelParameters     = modelParameters,
  modelError          = modelError,
  optimizer           = "FedorovWynnAlgorithm",
  optimizerParameters = list(
    elementaryProtocols   = initialElementaryProtocols,
    numberOfSubjects      = numberOfSubjects,
    proportionsOfSubjects = proportionsOfSubjects,
    showProcess           = FALSE
  ),
  designs             = list(designConstraint),
  fimType             = "population",
  outputs             = list("RespPK" = "Cc", "RespPD" = "E"),
  odeSolverParameters = list(atol = 1e-8, rtol = 1e-8)
)

optimizationFWPopFIM = run(optimizationFWPopFIM)
saveRDS(optimizationFWPopFIM,
        file.path(paths$data, "vignette_1_optimization_FedorovWynn_populationFIM.RDS"))

Display and plot Fedorov-Wynn results

Comparing RSE and D-criterion between evaluationPop and optimizationFWPopFIM quantifies the information gain from optimization with 30 vs 36 subjects.

For FedorovWynnAlgorithm, plotFrequencies() returns a bar chart showing how the 30 subjects are distributed across the selected elementary protocols. The number of non-zero bars is the support size of the D-optimal design.


show(optimizationFWPopFIM)
writeLines(capture.output(show(optimizationFWPopFIM)),
           file.path(paths$outputs, "vignette1_optimization_FedorovWynn_populationFIM_show.txt"))
fisherMatrix = getFisherMatrix(optimizationFWPopFIM)
getCorrelationMatrix(optimizationFWPopFIM)
getSE(optimizationFWPopFIM)
getRSE(optimizationFWPopFIM)
getShrinkage(optimizationFWPopFIM)
getDeterminant(optimizationFWPopFIM)
getDcriterion(optimizationFWPopFIM)

plotFWFrequencies = PFIM::plotFrequencies(optimizationFWPopFIM)
plotFW_SE = PFIM::plotSE(optimizationFWPopFIM)
plotFW_RSE = PFIM::plotRSE(optimizationFWPopFIM)
plotFWFrequencies
--- Optimal design ---

  Arms name Number of subjects Outcome  Dose     Sampling times
1      Arm1              12.06  RespPK    20    (0.25, 2, 4, 6)
2      Arm1              12.06  RespPD     .    (0.75, 2, 3, 6)
3      Arm2               7.37  RespPK 11.24 (0.25, 0.75, 1, 4)
4      Arm2               7.37  RespPD     .    (0.75, 2, 3, 6)
5      Arm3              10.57  RespPK    20 (0.25, 0.75, 1, 4)
6      Arm3              10.57  RespPD     .   (0.75, 2, 6, 12)

*************************************** 
  Population Fisher Matrix 
*************************************** 

                         μ_V          μ_Cl        μ_kout      μ_Imax      μ_IC50     μ_gamma         ω²_V        ω²_Cl      ω²_kout      ω²_Imax      ω²_IC50    ω²_gamma σ_slope_RespPK σ_inter_RespPD
μ_V             4.399120e+02  -32.44476547  6.479546e-05  -5.8026499  0.82567391 -0.54115176 0.000000e+00 0.000000e+00 0.000000e+00   0.00000000   0.00000000 0.000000000   0.000000e+00     0.00000000
μ_Cl           -3.244477e+01 1719.09514504 -6.305551e-02  -5.0522273  1.90909437  0.36447587 0.000000e+00 0.000000e+00 0.000000e+00   0.00000000   0.00000000 0.000000000   0.000000e+00     0.00000000
μ_kout          6.479546e-05   -0.06305551  8.724435e-01   0.6719275 -0.02841359 -0.04188020 0.000000e+00 0.000000e+00 0.000000e+00   0.00000000   0.00000000 0.000000000   0.000000e+00     0.00000000
μ_Imax         -5.802650e+00   -5.05222727  6.719275e-01 151.8820968 -3.58317296  4.35576466 0.000000e+00 0.000000e+00 0.000000e+00   0.00000000   0.00000000 0.000000000   0.000000e+00     0.00000000
μ_IC50          8.256739e-01    1.90909437 -2.841359e-02  -3.5831730  1.00784891 -0.02922741 0.000000e+00 0.000000e+00 0.000000e+00   0.00000000   0.00000000 0.000000000   0.000000e+00     0.00000000
μ_gamma        -5.411518e-01    0.36447587 -4.188020e-02   4.3557647 -0.02922741  0.96023665 0.000000e+00 0.000000e+00 0.000000e+00   0.00000000   0.00000000 0.000000000   0.000000e+00     0.00000000
ω²_V            0.000000e+00    0.00000000  0.000000e+00   0.0000000  0.00000000  0.00000000 9.685893e+02 1.265917e+00 3.143606e-05   0.24669467   0.76586183 0.020576534   2.204464e+02     0.01051979
ω²_Cl           0.000000e+00    0.00000000  0.000000e+00   0.0000000  0.00000000  0.00000000 1.265917e+00 3.028580e+02 2.113853e-04   0.02614577   0.41369145 0.003228195   3.927722e+01     0.01316118
ω²_kout         0.000000e+00    0.00000000  0.000000e+00   0.0000000  0.00000000  0.00000000 3.143606e-05 2.113853e-04 1.803158e+01   0.22420327   0.04847044 0.011481985   7.984847e-03     0.03852134
ω²_Imax         0.000000e+00    0.00000000  0.000000e+00   0.0000000  0.00000000  0.00000000 2.466947e-01 2.614577e-02 2.242033e-01 134.45077508  16.74357104 1.531084575   1.440512e+00     1.95850267
ω²_IC50         0.000000e+00    0.00000000  0.000000e+00   0.0000000  0.00000000  0.00000000 7.658618e-01 4.136914e-01 4.847044e-02  16.74357104 126.35481896 0.093799691   1.226609e+01     2.50237286
ω²_gamma        0.000000e+00    0.00000000  0.000000e+00   0.0000000  0.00000000  0.00000000 2.057653e-02 3.228195e-03 1.148199e-02   1.53108457   0.09379969 0.914300103   1.741369e-01     0.10479196
σ_slope_RespPK  0.000000e+00    0.00000000  0.000000e+00   0.0000000  0.00000000  0.00000000 2.204464e+02 3.927722e+01 7.984847e-03   1.44051173  12.26608742 0.174136864   2.800700e+03     0.34850788
σ_inter_RespPD  0.000000e+00    0.00000000  0.000000e+00   0.0000000  0.00000000  0.00000000 1.051979e-02 1.316118e-02 3.852134e-02   1.95850267   2.50237286 0.104791955   3.485079e-01     0.42642687

*************************************** 
  Fixed effects (μ) 
*************************************** 

                  μ_V          μ_Cl        μ_kout      μ_Imax      μ_IC50     μ_gamma
μ_V      4.399120e+02  -32.44476547  6.479546e-05  -5.8026499  0.82567391 -0.54115176
μ_Cl    -3.244477e+01 1719.09514504 -6.305551e-02  -5.0522273  1.90909437  0.36447587
μ_kout   6.479546e-05   -0.06305551  8.724435e-01   0.6719275 -0.02841359 -0.04188020
μ_Imax  -5.802650e+00   -5.05222727  6.719275e-01 151.8820968 -3.58317296  4.35576466
μ_IC50   8.256739e-01    1.90909437 -2.841359e-02  -3.5831730  1.00784891 -0.02922741
μ_gamma -5.411518e-01    0.36447587 -4.188020e-02   4.3557647 -0.02922741  0.96023665

*************************************** 
  Variance components (ω², γ², σ) 
*************************************** 

                       ω²_V        ω²_Cl      ω²_kout      ω²_Imax      ω²_IC50    ω²_gamma σ_slope_RespPK σ_inter_RespPD
ω²_V           9.685893e+02 1.265917e+00 3.143606e-05   0.24669467   0.76586183 0.020576534   2.204464e+02     0.01051979
ω²_Cl          1.265917e+00 3.028580e+02 2.113853e-04   0.02614577   0.41369145 0.003228195   3.927722e+01     0.01316118
ω²_kout        3.143606e-05 2.113853e-04 1.803158e+01   0.22420327   0.04847044 0.011481985   7.984847e-03     0.03852134
ω²_Imax        2.466947e-01 2.614577e-02 2.242033e-01 134.45077508  16.74357104 1.531084575   1.440512e+00     1.95850267
ω²_IC50        7.658618e-01 4.136914e-01 4.847044e-02  16.74357104 126.35481896 0.093799691   1.226609e+01     2.50237286
ω²_gamma       2.057653e-02 3.228195e-03 1.148199e-02   1.53108457   0.09379969 0.914300103   1.741369e-01     0.10479196
σ_slope_RespPK 2.204464e+02 3.927722e+01 7.984847e-03   1.44051173  12.26608742 0.174136864   2.800700e+03     0.34850788
σ_inter_RespPD 1.051979e-02 1.316118e-02 3.852134e-02   1.95850267   2.50237286 0.104791955   3.485079e-01     0.42642687

********************************************* 
  Determinant, condition numbers and D-criterion 
 *********************************************** 

Determinant: 5.739199e+21 
D-criterion: 35.82644 
Condition number (fixed effects): 2239.911 
Condition number (variance components): 8224.375 

*************************************** 
  Parameters estimation 
*************************************** 

Parameter               Value           SE     RSE(%)
μ_V                 0.740000   0.04776611   6.454880
μ_Cl                0.280000   0.02416351   8.629825
μ_kout              6.140000   1.07524517  17.512136
μ_Imax              0.760000   0.09142652  12.029805
μ_IC50              9.220000   1.04615256  11.346557
μ_gamma             2.770000   1.10108306  39.750291
ω²_V               0.099856   0.03242338  32.470135
ω²_Cl              0.207936   0.05751467  27.659796
ω²_kout            0.896809   0.23552053  26.262062
ω²_Imax            0.192721   0.08983672  46.614912
ω²_IC50            0.204304   0.09489532  46.448095
ω²_gamma           3.101121   1.06789153  34.435662
σ_slope_RespPK       0.210000   0.01908902   9.090012
σ_inter_RespPD       9.600000   1.69304377  17.635873

[1] 35.82644


outputFile = "Example01_OptimizationFWPopFIM.html"
Report(optimizationFWPopFIM, paths$reports, outputFile, plotOptions)

Multiplicative algorithm

The Multiplicative algorithm (cocktail / multiplicative weights update) is an iterative continuous relaxation approach to D-optimal design in a discrete candidate space. It proceeds as:

  1. Initialize uniform weights wâ‚– = 1/K for each of K candidate protocols.
  2. At each iteration, update weights multiplicatively: wₖ(t+1) ∝ wₖ(t) × dₖ(ξₜ), where dₖ(ξₜ) is the normalized directional derivative of the D-criterion for protocol k under the current design ξₜ.
  3. After convergence (change in D-criterion < delta), set to zero all weights below weightThreshold, yielding a sparse approximate D-optimal design.

Unlike the Fedorov-Wynn algorithm, the MA explores the full candidate space simultaneously without requiring initial elementary protocols, and may converge to a different local optimum.

MultiplicativeAlgorithm optimizerParameters:

Parameter Description
lambda step-size dampening factor (0 < λ < 1); values close to 1 slow convergence but reduce oscillations
numberOfIterations maximum number of multiplicative update cycles
weightThreshold protocols with weight < threshold at convergence are zeroed (here 0.01 = < 1% of subjects)
delta convergence tolerance on the relative D-criterion change (here 1e-4: stop when improvement < 0.01%)
showProcess logical; if TRUE, prints per-iteration D-criterion values 

optimizationMultPopFIM = Optimization(
  name                = "PKPD_ODE_multi_doses_populationFIM",
  modelEquations      = modelEquations,
  modelParameters     = modelParameters,
  modelError          = modelError,
  optimizer           = "MultiplicativeAlgorithm",
  optimizerParameters = list(
    lambda             = 0.99,
    numberOfIterations = 1000,
    weightThreshold    = 0.01,
    delta              = 1e-04,
    showProcess        = FALSE
  ),
  designs             = list(designConstraint),
  fimType             = "population",
  outputs             = list("RespPK" = "Cc", "RespPD" = "E"),
  odeSolverParameters = list(atol = 1e-8, rtol = 1e-8)
)

optimizationMultPopFIM = run(optimizationMultPopFIM)
saveRDS(optimizationMultPopFIM,
        file.path(paths$data, "vignette_1_optimization_multiplicativeAlgorithm_populationFIM.RDS"))

Display and plot Multiplicative algorithm results

For MultiplicativeAlgorithm, plotWeights() returns a bar chart of the final weight of each candidate protocol after convergence. Non-zero weights define the design support; comparing this plot to plotFrequencies() from FW reveals whether both algorithms converge to the same support, confirming robustness.


show(optimizationMultPopFIM)
writeLines(capture.output(show(optimizationMultPopFIM)),
           file.path(paths$outputs, "vignette1_optimization_MultiplicativeAlgorithm_populationFIM_show.txt"))
fisherMatrix = getFisherMatrix(optimizationMultPopFIM)
getCorrelationMatrix(optimizationMultPopFIM)
getSE(optimizationMultPopFIM)
getRSE(optimizationMultPopFIM)
getShrinkage(optimizationMultPopFIM)
getDeterminant(optimizationMultPopFIM)
getDcriterion(optimizationMultPopFIM)

plotMultWeights = PFIM::plotWeights(optimizationMultPopFIM)
plotMult_SE = PFIM::plotSE(optimizationMultPopFIM)
plotMult_RSE = PFIM::plotRSE(optimizationMultPopFIM)
ggsave(file.path(paths$figures, "vignette1_optimization_MultiplicativeAlgorithm_populationFIM_weights.pdf"),
       plotMultWeights, width = 8, height = 5)
plotMultWeights
--- Optimal design ---

  Arms name Number of subjects Outcome  Dose     Sampling times
1    Arm456              10.81  RespPK    20    (0.25, 2, 4, 6)
2    Arm456              10.81  RespPD     .    (0.75, 2, 3, 6)
3    Arm462              10.31  RespPK    20 (0.25, 0.75, 1, 4)
4    Arm462              10.31  RespPD     .   (0.75, 2, 6, 12)
5    Arm368               7.39  RespPK 11.24 (0.25, 0.75, 1, 4)
6    Arm368               7.39  RespPD     .    (0.75, 2, 3, 6)
7    Arm451               1.49  RespPK    20 (0.25, 0.75, 1, 4)
8    Arm451               1.49  RespPD     .    (0.75, 2, 3, 6)

--- Optimal mixture weights ---

 Grid cell Weight N subjects
       368 0.3602         11
       451 0.3437         10
       456 0.2465          7
       462 0.0496          2

  (N subjects = Hamilton / largest-remainder of 30 * weight; sum(N) = 30)

*************************************** 
  Population Fisher Matrix 
*************************************** 

                        μ_V          μ_Cl      μ_kout      μ_Imax      μ_IC50     μ_gamma         ω²_V        ω²_Cl      ω²_kout      ω²_Imax      ω²_IC50    ω²_gamma σ_slope_RespPK σ_inter_RespPD
μ_V            441.31521821  -30.19202429  0.00216574  -5.7727396  0.77506468 -0.54061196 0.000000e+00 0.000000e+00 0.000000e+00   0.00000000   0.00000000 0.000000000   0.000000e+00    0.000000000
μ_Cl           -30.19202429 1716.54061463 -0.06801023  -4.9118108  1.98401441  0.38757590 0.000000e+00 0.000000e+00 0.000000e+00   0.00000000   0.00000000 0.000000000   0.000000e+00    0.000000000
μ_kout           0.00216574   -0.06801023  0.87229393   0.6788107 -0.02842157 -0.04223823 0.000000e+00 0.000000e+00 0.000000e+00   0.00000000   0.00000000 0.000000000   0.000000e+00    0.000000000
μ_Imax          -5.77273961   -4.91181077  0.67881072 151.5036810 -3.58077140  4.37444810 0.000000e+00 0.000000e+00 0.000000e+00   0.00000000   0.00000000 0.000000000   0.000000e+00    0.000000000
μ_IC50           0.77506468    1.98401441 -0.02842157  -3.5807714  1.00932203 -0.02925790 0.000000e+00 0.000000e+00 0.000000e+00   0.00000000   0.00000000 0.000000000   0.000000e+00    0.000000000
μ_gamma         -0.54061196    0.38757590 -0.04223823   4.3744481 -0.02925790  0.95933219 0.000000e+00 0.000000e+00 0.000000e+00   0.00000000   0.00000000 0.000000000   0.000000e+00    0.000000000
ω²_V             0.00000000    0.00000000  0.00000000   0.0000000  0.00000000  0.00000000 9.747086e+02 1.143241e+00 5.904948e-05   0.24533182   0.71143518 0.020534866   2.186688e+02    0.009852846
ω²_Cl            0.00000000    0.00000000  0.00000000   0.0000000  0.00000000  0.00000000 1.143241e+00 3.019533e+02 2.670295e-04   0.02588491   0.45771094 0.003320576   4.011414e+01    0.013799621
ω²_kout          0.00000000    0.00000000  0.00000000   0.0000000  0.00000000  0.00000000 5.904948e-05 2.670295e-04 1.802540e+01   0.22796734   0.04839936 0.011640400   9.095069e-03    0.038809606
ω²_Imax          0.00000000    0.00000000  0.00000000   0.0000000  0.00000000  0.00000000 2.453318e-01 2.588491e-02 2.279673e-01 133.75270473  16.75419698 1.544249521   1.440902e+00    1.954132465
ω²_IC50          0.00000000    0.00000000  0.00000000   0.0000000  0.00000000  0.00000000 7.114352e-01 4.577109e-01 4.839936e-02  16.75419698 126.61505699 0.093824817   1.297849e+01    2.487552632
ω²_gamma         0.00000000    0.00000000  0.00000000   0.0000000  0.00000000  0.00000000 2.053487e-02 3.320576e-03 1.164040e-02   1.54424952   0.09382482 0.912709894   1.767293e-01    0.104477346
σ_slope_RespPK   0.00000000    0.00000000  0.00000000   0.0000000  0.00000000  0.00000000 2.186688e+02 4.011414e+01 9.095069e-03   1.44090233  12.97849240 0.176729295   2.799473e+03    0.358466134
σ_inter_RespPD   0.00000000    0.00000000  0.00000000   0.0000000  0.00000000  0.00000000 9.852846e-03 1.379962e-02 3.880961e-02   1.95413247   2.48755263 0.104477346   3.584661e-01    0.427005132

*************************************** 
  Fixed effects (μ) 
*************************************** 

                 μ_V          μ_Cl      μ_kout      μ_Imax      μ_IC50     μ_gamma
μ_V     441.31521821  -30.19202429  0.00216574  -5.7727396  0.77506468 -0.54061196
μ_Cl    -30.19202429 1716.54061463 -0.06801023  -4.9118108  1.98401441  0.38757590
μ_kout    0.00216574   -0.06801023  0.87229393   0.6788107 -0.02842157 -0.04223823
μ_Imax   -5.77273961   -4.91181077  0.67881072 151.5036810 -3.58077140  4.37444810
μ_IC50    0.77506468    1.98401441 -0.02842157  -3.5807714  1.00932203 -0.02925790
μ_gamma  -0.54061196    0.38757590 -0.04223823   4.3744481 -0.02925790  0.95933219

*************************************** 
  Variance components (ω², γ², σ) 
*************************************** 

                       ω²_V        ω²_Cl      ω²_kout      ω²_Imax      ω²_IC50    ω²_gamma σ_slope_RespPK σ_inter_RespPD
ω²_V           9.747086e+02 1.143241e+00 5.904948e-05   0.24533182   0.71143518 0.020534866   2.186688e+02    0.009852846
ω²_Cl          1.143241e+00 3.019533e+02 2.670295e-04   0.02588491   0.45771094 0.003320576   4.011414e+01    0.013799621
ω²_kout        5.904948e-05 2.670295e-04 1.802540e+01   0.22796734   0.04839936 0.011640400   9.095069e-03    0.038809606
ω²_Imax        2.453318e-01 2.588491e-02 2.279673e-01 133.75270473  16.75419698 1.544249521   1.440902e+00    1.954132465
ω²_IC50        7.114352e-01 4.577109e-01 4.839936e-02  16.75419698 126.61505699 0.093824817   1.297849e+01    2.487552632
ω²_gamma       2.053487e-02 3.320576e-03 1.164040e-02   1.54424952   0.09382482 0.912709894   1.767293e-01    0.104477346
σ_slope_RespPK 2.186688e+02 4.011414e+01 9.095069e-03   1.44090233  12.97849240 0.176729295   2.799473e+03    0.358466134
σ_inter_RespPD 9.852846e-03 1.379962e-02 3.880961e-02   1.95413247   2.48755263 0.104477346   3.584661e-01    0.427005132

********************************************* 
  Determinant, condition numbers and D-criterion 
 *********************************************** 

Determinant: 5.731602e+21 
D-criterion: 35.82305 
Condition number (fixed effects): 2243.102 
Condition number (variance components): 8189.148 

*************************************** 
  Parameters estimation 
*************************************** 

Parameter               Value           SE     RSE(%)
μ_V                 0.740000   0.04768087   6.443361
μ_Cl                0.280000   0.02418124   8.636156
μ_kout              6.140000   1.07544879  17.515453
μ_Imax              0.760000   0.09163559  12.057314
μ_IC50              9.220000   1.04542683  11.338686
μ_gamma             2.770000   1.10273794  39.810034
ω²_V               0.099856   0.03231505  32.361653
ω²_Cl              0.207936   0.05760329  27.702412
ω²_kout            0.896809   0.23556120  26.266596
ω²_Imax            0.192721   0.09009104  46.746872
ω²_IC50            0.204304   0.09470714  46.355989
ω²_gamma           3.101121   1.06890883  34.468466
σ_slope_RespPK       0.210000   0.01909074   9.090826
σ_inter_RespPD       9.600000   1.68992464  17.603382

[1] 35.82305


outputFile = "Example01_OptimizationMultPopFIM.html"
Report(optimizationMultPopFIM, paths$reports, outputFile, plotOptions)

References

Flores-Murrieta, Francisco, Holly Kimko, Dora Flores-Acevedo, Francisco López-Muñoz, William Jusko, Mark Sale, and Gilberto Castañeda-Hernández. 1998. “Pharmacokinetic–Pharmacodynamic Modeling of Tolmetin Antinociceptive Effect in the Rat Using an Indirect Response Model: A Population Approach.” Journal of Pharmacokinetics and Biopharmaceutics 26 (November): 547–57.