Design Evaluation and Optimization of Population Pharmacokinetics Model Using an R Package PopED

Abstract

The rising cost in drug development has led to continuous calls for more efficient and powerful design optimization and analysis tools for pkpd studies. More sophisticated models are increasingly used to reflect reality, and current models invariably include non-linear mixed effects models that frequently require specialized computational tools for the design and analysis of the study. Population PK analysis employs non-linear mixed effects models and evaluates designs and optimizes them via a specialized yet versatile software package called PopED (version 0.6.0) in programming language R (version 4.2.2) for pharmacometrics analyses. We demonstrate the utilities of the package when different models and statistical criteria are used in real settings to determine the optimal sampling times and optimal dose levels for the subjects. We provide two applications; the first is illustrative and the second is a new application on developing an optimal dosing scheme for a two-compartment PK model with perturbation. Our target audiences are mathematicians and statisticians who are not aware of this useful and powerful analytic tool.

1. Introduction

Finding the best dose regimen and sampling times to ascertain the pharmacokinetics and pharmacodynamics of a drug are important aspects of drug development [1]. The cost of drug development continues to skyrocket, and design issues have been increasingly emphasized, not only for reining in cost but also for statistical efficiency and ethical issues for patients [2]. Finding the best design for the most accurate statistical inference at minimal cost is therefore a particularly pertinent issue in drug development.

Pharmacokinetics (PK) models delineate how the drug, once it is consumed, interacts with our physiological system, from the administration to the complete elimination of the drug from the body [3]. PK investigates how the substances in the drug change through the course of its being absorbed, distributed, metabolized, and excreted (often referred to as the ADME process from the first letters of the four verbs). PK often focuses on, but is not limited to, dose–concentration relationships, where we assume that there is a mathematical expression that describes how different dose levels affect drug concentration in the body.

The traditional individual approach to understanding PK typically focuses on a relatively homogeneous group of subjects with only a single factor most likely to affect the PK and pharmacodynamics (PD) of the drug. Broadly speaking, PD is the science on how the body reacts to drugs. Naturally, the PK relationships of a drug depend on the individual traits, such as the demographic, physiologic, and genetic characteristics, of each person. Weight, renal function, and genetic predispositions are some examples of covariates that affect the PK relationship. One drawback of traditional PK is that it cannot properly account for the variability of such factors, nor the variability arising among individuals.

One way to account for the variability among the individuals is to use the population approach, which has several merits over the traditional approach. The former typically requires denser and balanced data points from each individual to estimate the model parameters well, but the population PK approach does not have such requirements and have been shown to be generally more robust when we have unbalanced data. Second, the population PK approach accounts for one or more of the covariate effects, and third, the population PK approach has random effects, which can account for the unknown variability stemming from inherent individual differences that are not explained by having fixed covariates. The upshot is that the population approach is more flexible and provides a holistic analysis of the PK relationship. As such, it has been widely used in the pharmaceutical industry for drug development via clinical trials since the 1990s [4]. They provide critical clinical information of the drug regarding the optimal dosage regimen and population-specific covariate effects. This suggests that in practice, we should adopt the population approach and find efficient, if not optimal, designs for the PK models that provide maximal information for designing optimal dose regimens, sampling times, and estimating population specific covariate parameters at minimum cost [5].

In the clinical trials Phase I, Phase II, and Phase III, each phase requires more patients than the previous phase. Across these trial phases, we seek efficient, if not optimal, designs for the PK models, which is equivalent to finding a design that provides maximal information for designing optimal dose regimens, sampling times, and estimating population-specific covariate parameters at minimum cost [5]. Without the optimized doses and sampling plans that maximize the informativeness of the experiments, recruiting participants for dense sampling is both costly and unnecessary. Therefore, organizing and evaluating designs for PK analysis based on good evaluation and optimization methods are crucial components in drug development.

For this reason, it is necessary to have an efficient and user-friendly computational tool for design evaluation and optimization. Some nascent recognition of importance of such a step is apparent in the Food and Drug Administration’s suggestion, through its 1999 Guidance for Industry Population Pharmacokinetics, that the design evaluation simulations be performed for possible clinical trial designs [6]. Historically, it is common to assume that a statistical model is given and an optimal design is found for a particular objective under the model assumptions. In recent years, oncologists and medical researchers have also embraced model-based designs [7].

Following convention, the worth of a design is measured by its Fisher information matrix (FIM), which is the expectation of the second derivatives of the total log-likelihood function with respect to the model parameters, apart from a multiplicative constant. This matrix depends on the model, the design that generated the data and the unknown parameters of the nonlinear model. A common design criterion is to find a design that maximizes the determinant of the information matrix among all designs in the given compact interval, which may be the sampling time range or the dose range. Because of the presence of the unknown parameters in the matrix which we wish to estimate, we cannot optimize the determinant directly. Nominal values for the unknown parameters are best-guess values for the unknown parameters, and they may come from experts’ opinions or previous similar trails. After substituting the nominal values into the information matrix, the objective function (i.e., the determinant) depends only the design, and we can then optimize it by finding the optimal number of design points where they are and the proportion of observations to take at each of the points. The resulting optimal design is local because it depends on the nominal values for the unknown parameters. The optimization itself can be difficult when there are many interacting factors in the model or there are random components in the model, along with constraints on the longitudinal study. For example, requiring the number and locations of the sampling time points can become very complicated quickly, and the information matrix is computed numerically before it is optimized according to some design criteria. Kiefer (1980) formulated all objectives as convex functions of the information matrix so that the optimization problem becomes convex, and tools from the convex analysis can be applied.

Software for performing such a task started in Splus, and MATLAB including self-created software like those from [8,9,10]. Other software tools, such as PFIM, PopDes, or PopED [11,12], are also actively used in population PK design evaluation and optimization [6], with the more recent one being the dose-optimizing package OptiDose from the nonlinear mixed effects modeling software NONMEM [13,14]. This is in alignment with growing interest in utilizing statistical optimization methodologies in earlier stages of drug development [15]. In the realm of population PK design, FIM-based evaluation and optimization is well studied [16,17,18,19], including traditional methodologies using other optimization criteria for population PK [20,21]. With growing interest in the programming language R, there is also a need for having a R-based tool for the design evaluation and optimization in population PK.

The purpose of this paper is to present a useful computing package called PopED for statisticians analyzing NLMEM population PK models. PopED uses the programming language R for the design optimization of non-linear mixed effects PK models, which has the potential to be widely used in the statistical community. This article demonstrates how the PopED package is used in PK design evaluation and optimization via two applications: (i) a one-compartment first-order PK model with single oral dosing and multiple bolus administration, and (ii) a single dosing two-compartment PK model with perturbation. The optimization simulation will follow the mathematical exposition of the PK models for each example to provide the readers with a more comprehensive understanding of the models. To this end, we first give a brief introduction of non-linear mixed effects models with some technical details for the models in the two applications. We then use PopED to design and evaluate various designs via the simulation tool in the software. The paper concludes with a discussion. The simulation code can be found in https://github.com/HowonRyu/PK_PopED.

2. Statistical Models

This section first reviews non-linear mixed effects models (NLMEMs) and their structure with two applications for PK modeling. The first application uses a one-compartment first-order absorption PK model, and the second application uses a single dosing two-compartment pharmacokinetics model with perturbation. The first is illustrative, and the second is a new application. We then describe how to find optimal designs and how to evaluate whether a design is optimal or not by ascertaining the design efficiency measures.

2.1. Non-Linear Mixed Effects Models for Population Pharmacokinetics

NLMEMs are commonly used in PK studies because (i) PK models are invariably nonlinear, and (ii) population PK studies have repeated measures data, where each individual is treated with a dose of a drug and observed multiple times over the specified study period [22]. Further, the random effects in the model add an extra layer of variability into the population PK models.

NLMEMs, as the name suggests, have both fixed and random effects in the models. Fixed effects are usually a subject-specific characteristic that does not change over the study period where repeated measures are observed. Random effects are variables assigned to different unit groups of the repeated measures to account for the differences in between the groups [22]. Random effects usually account for unspecified variability stemming from differences that are not explained by the differences in fixed-effect covariates.

The setup for NLMEMs assumes that there are N individuals for the study, where N is predetermined, and the i-th subject has $n_i$ repeated measures $Y_{ij}$ at time point $t_{ij},j=1,\dots ,n_i$ and $i=1,\dots ,N$. Individual factors, such as age, weight, or blood pressure, are integrated into the modeling equation as covariates, although NLMEM PK analysis is possible without these covariates. The model structure varies depending on the application, and we assume that data from different subjects are independent. A general form of such a model is as follows:

• The observed data for the i-th individual: ${(Y_i,U_i,A_i)}^{\prime }$ for $i=1,\dots ,N.$
• $Y_i={(Y_{i1},\dots ,Y_{i{n}_i})}^{\prime },$ $U_i=(U_{i1},\dots ,$$U_{i{n}_i}{)}^{\prime },$ and $A_i={(A_i,\dots ,A_N)}^{\prime }.$
• $Y_{ij}=m(t_{ij},U_i,{\theta }_i)+{ϵ}_{ij}$, where ${ϵ}_{ij}\sim N(0,{\sigma }^2)$ and $j=1,\dots ,n_i.$
• ${\theta }_i=d(A_i,\beta ,{\eta }_i)$ where ${\eta }_{\sim }MVN(0,\mathbf{G})$ and the ${\eta }_i$’s are uncorrelated.

Here and throughout, each outcome $Y_{ij}$ is modeled with three factors: the time component $t_{ij}$, the within-individual (intra-individual) covariate $U_i$, and the between-individual (inter-individual) component ${\theta }_i$. An example of $U_i$ as a covariate is the drug dose, which can have multiple values for an individual. Examples of ${\theta }_i$ as covariates are the absorption rate constant ($K_a$), clearance ($CL$), or volume of distribution (V), which will be further explained in this section. The vector $A_i$ usually contains subject-character covariates, such as weight, blood pressure, or age. The functions m and d have different forms depending on the PK model. The parameter ${\theta }_i$ typically comprises three factors: $A_i$, the fixed effects; $\beta$, the fixed-effect coefficients; and ${\eta }_i$, the random effects. The residual error ${ϵ}_{ij}$ accounts for the residual variability for $Y_{ij}$, and it is assumed to be normally distributed with mean 0 and variance ${\sigma }^2$. The random-effect $\eta$ follows a multivariate normal distribution with mean 0 and covariance matrix $\mathbf{G}$, and they are often assumed to be uncorrelated with one another, implying that their covariance matrix is simply a diagonal matrix $diag(var\left({\eta }_1\right),\dots ,var\left({\eta }_m\right))$ when assuming m different parameters. The variables $Y_{ij}$, $A_i$, $U_i$, and $t_{ij}$ are assumed to be known or observable. The unknown parameters to be estimated are $\beta$, ${\eta }_i$, ${ϵ}_{ij}$, ${\sigma }^2$, and $\mathbf{G}$.

We next use two PK models to demonstrate how NLMEMs are formulated and used to analyze data from a longitudinal study using PopED.

2.2. Application 1: One-Compartment First-Order Absorption Single and Multiple Dosing PK Model

The first application is taken from the instructional webpage of PopED created by Andrew C. Hooker (https://andrewhooker.github.io/PopED/). This application is a one-compartment first-order absorption PK model with single oral dosing and multiple IV bolus administration (intravenous injection). A compartment is a conceptual unit within the human body, where the organs and tissues belonging to that compartment are thought to have similar drug distribution. Therefore, the one-compartment model assumes instant distribution of the drug upon administration, that the drug distribution is similar throughout the body, and that there is no secondary distribution of the drug in peripheral parts of the main compartment [23]. In a similar way, the two-compartment model has two such systems, with a likely slower distribution of the drug. The first-order absorption model is also called a linear absorption model because the model assumes a constant absorption rate. It supposes that the relationship between the dose level and the subsequent change of drug concentration is linear. The mathematical model is

$$\begin{array}{cc}\hfill Y=\frac{D·F}{V}\frac{Ka}{Ka-\frac{CL}{V}}·& \{exp(\frac{-CL}{V}(t-(N_d-1)\tau ))\frac{1-exp(-N_d\frac{CL}{V}\tau )}{1-exp(-\frac{CL}{V}\tau )}\hfill \\ \hfill & -exp(-K_{a}t-(N_d-1)\tau )\frac{1-exp(-N_d{K}_a\tau )}{1-exp(-K_a\tau )}\}.\hfill \end{array}$$

(1)

In (1), the outcome variable Y (mmol/L) is the drug plasma concentration, and the variable t (h) denotes the time. D (mmol) denotes dose, and there is parameter $\tau$ (h) which denotes dose interval since Application 1 is a multiple dosing model. The variable F stands for bioavailability, which is the fraction of the drug absorbed by the systemic circulation, and V (L) is the volume of distribution of the drug. The variable $K_a$ represents an absorption rate constant, and $CL$ (L/h) is clearance, which is the volume of the plasma that is cleared of the drug. Finally, we let $N^{\ast }$ be the largest integer that is not greater than the value calculated by diving t (time) by $\tau$ (dose interval), and let $N_d=N^{\ast }+1$.

This PK model can be formulated using a general form of a PK model given by the following:

• $Y_{ij}=m(t_{ij},U_i,{\theta }_i)+{ϵ}_{ij},$ where ${ϵ}_{ij}\sim N(0,{\sigma }^2),j=1,\dots ,n_i$, and m is presented as (1).
• ${\theta }_i=(V_i,K_{ai},C{L}_i,F_i)=d(A_i,\beta ,{\eta }_i)$ for $i=1,\dots ,N,$ where
  $V_i={\beta }_{0V}·e^{{\eta }_{i1}}=\overline{V}·e^{{\eta }_{i1}}$, $K_{ai}={\beta }_{0K_a}·e^{{\eta }_{i2}}=\overline{K_a}·e^{{\eta }_{i2}}$, $C{L}_i={\beta }_{0_{CL}}·e^{{\eta }_{i3}}=\overline{CL}·e^{{\eta }_{i3}}$, and $F_i={\beta }_0=\overline{F}$.
• $U_{ij}={(D_{ij},{\tau }_{ij})}^{\prime },i=1,\dots ;N,j=1,\dots ,n_i.$

The covariates ${\theta }_i$ are the volume of distribution, absorption rate constant, clearance, and bioavailability. This model does not have any subject characteristic covariate ($A_i$), and its fixed effects are linearly modeled. Therefore, ${\theta }_i=(V_i,K_{a_i},C{L}_i,F_i)$ is a combination of mean fixed effects and the random effects $({\eta }_{i1},{\eta }_{i2},{\eta }_{i3})$. This is a special case of the more general expression $V_i=({\beta }_{0V}+{\beta }_{1V}·A_i)·e^{{\eta }_{i1}}$, where $A_i=0$, and the intercept ${\beta }_{0V}$ is replaced with $\overline{V}$, which stands for the averaged volume across all patients. The same holds for $K_{a_i}$ and $C{L}_i$. We note that $F_i$ is considered fixed and has no random-effect component. All the random effects follow an exponential model with a lognormal distribution.

2.3. Application 2: Single Dosing Two-Compartment PK Model with Perturbation

This application concerns single dosing two-compartment PK models with design optimization for perturbation PK models [24]. Unlike the first application, such a model involves a hemodialysis as perturbation during the dosing cycle as the intermittent exogenous elimination factor of the drug from blood. The model takes into account the loss of drug concentration due to perturbation, and we denote the endogenous elimination rate and the perturbation rate of the drug by $k_{10}$ and $k_D$, respectively. This application supposes a two-compartment model, where there is a distribution of the drug between central and peripheral compartment, which is another factor that influences the PK dynamics.

Specifically, the model assumes single administration of a drug by intravenous infusion to the central compartment ($k_R$). The dialytic clearance ($k_D$) happens mid-cycle, and the regular endogenous elimination of the drug ($k_{10}$) is assumed for the remaining cycle. In the above equations, $c_1$ and $c_2$ are the drug concentrations in the central and peripheral compartments, respectively, and $v_1$ is the central compartment volume. The parameters $k_{12}$ and $k_{21}$ are the rates at which the drug is distributed between the central and peripheral compartments (L/h). The distribution rate from the central to the peripheral compartments is $k_{12}$, and the distribution rate from the peripheral to the central compartment is $k_{21}$. Mathematically, the model is defined by the two differential equations, which were proposed and solved analytically by Shotwell (2016) [24]:

$$\frac{d{c}_1}{dt}=-k_{12}{c}_1+k_{22}{c}_2-k_{10}{c}_1+\frac{k_R}{v_1}-\frac{k_D}{v_1}{c}_1,$$

(2)

$$\frac{d{c}_2}{dt}=k_{12}{c}_1-k_{21}{c}_2.$$

(3)

Shotwell (2016) presents a piecewise analytic solution, where the drug concentration $Y_{ij}$ is represented by a function of t, $k_{10}$, $k_{12}$, $k_{21}$, and $v_1$, which are the main between-individual parameters of the model. The parameters $k_D$, $c_1$ and $c_2$ are not explicit in the analytic solution equation of Shotwell (2016) [24], because $k_D$ is a function of time t and $c_1$, while $c_2$ is of the rest of the parameters. The parameter $k_R$ is the dose of the drug infusion, which is usually predetermined, and so it is considered a constant. The model used by Shotwell may be restructured as follows:

• $Y_{ij}=m(t_{ij},U_i,{\theta }_i)+{ϵ}_{ij},$ where ${ϵ}_{ij}\sim N(0,{\sigma }^2),j=1,\dots ,n_i,$ and m is the model suggested by the analytical solution obtained from (2) and (3).
• ${\theta }_i=d(A_i,\beta ,{\eta }_i)=(k_{10i},k_{12i},k_{21i},v_{1i})$ for $i=1,\dots ,N,$ where $k_{10i}={\beta }_{0k_{10}}·e^{{\eta }_{i1}}={\tilde{k}}_{10}·e^{{\eta }_{i1}}$, $k_{12i}={\beta }_{0k_{12}}·e^{{\eta }_{i2}}={\tilde{k}}_{12}·e^{{\eta }_{i2}}$, $k_{21i}={\beta }_{0k_{21}}·e^{{\eta }_{i3}}={\tilde{k}}_{21}·e^{{\eta }_{i3}}$, and $v_{1i}={\beta }_{0v_1}·e^{{\eta }_{i4}}={\tilde{v}}_1·e^{{\eta }_{i4}}$.
• $U_{ij}=k_R,$ for $i=1,\dots ,N,$ and $j=1,\dots ,n_i.$

The between-individual covariates ${\theta }_i$ comprise $(k_{10i},k_{12i},k_{21i},v_{1i})$ and they represent, respectively, the endogenous clearance, distribution rate from the central to the peripheral compartment, distribution rate from the peripheral to the central compartment, and the central compartment volume. As is the case in Application 1, this model does not have any subject covariate ($A_i$). Because the covariates are linearly modeled, each component is defined as a combination of population median fixed effects $({\tilde{k}}_{10i},{\tilde{k}}_{12i},{\tilde{k}}_{21i},{\tilde{v}}_{1i})$, and the random effects $({\eta }_{i1},{\eta }_{i2},{\eta }_{i3},{\eta }_{i4})$. The above expression is obtained from the general expression $k_{10i}=({\beta }_{0k_{10}}+{\beta }_{1k_{10}}·A_i)·e^{{\eta }_{i1}}$ when $A_i=0$ after the intercept ${\beta }_{0k_{10}}$ is replaced by ${\tilde{k}}_{10}$, the median endogenous clearance rate of the population. The same holds for the rest of ${\theta }_i$. All random effects follow an exponential model with a lognormal distribution.

3. Design Optimization Method

This section discusses the fundamentals of constructing optimal designs, how to find them and confirm that they are indeed optimal. To fix ideas, we work with a relatively simple model to illustrate the concepts and difficulty involved.

3.1. Design Optimization

To appreciate the utility and capabilities of the software PopED, it is instructive to review the basics of finding an optimal design and how to confirm a design is optimal in an idealized case when we have a large sample and there are no random effects in the model. For simplicity, assume we have a linear regression model defined over a compact design space S and errors are independent, normally and identically distributed with mean 0 and constant variance. The design space for PK/PD studies is usually the range of sampling times or the range of dose concentrations allowed in the study. A usual goal is to find a design that estimates the vector of coefficients in the mean function as accurately as possible.

Designs have two forms: approximate and exact. Approximate designs are probability measures defined on S, and they are characterized by the number of design points (which are sampling time points, or dose levels), their locations in S and the proportion ($p_i$) of observations to take at each of these points. If the study is allowed to have a total of N observations, then we first round up each $N\ast p_i$ to an integer and take that number of observations at the i-th location. Exact designs assume that the total number of observations is pre-specified, and the design problem is similar, except the number of replicates at the locations are directly optimized.

Following convention, we measure the worth of a design $\xi$ by its information matrix $M(\xi ,\theta )$, which is proportional to the expectation of the second derivatives of the total log-likelihood function with respect to the parameters. Here, $\theta$ is the unknown parameter in the function, which enters into the information matrix when the model is nonlinear. The design criterion or objective function is $\varphi (\xi ,\theta )=-log\left|M\right(\xi ,\theta \left)\right|$, and we want to find a design $\xi$ (exact or approximate) to optimize the objective function over all designs on S. It can be shown that if the focus is on approximate designs, for fixed $\theta$, $\varphi$ is a convex function on the set of all $\xi$’s on S and consequently, the design problem is a convex optimization problem. Results from the convex analysis can then be usefully employed to find the optimal approximate designs and assess a design from the optimum, without knowing the latter, using the directional derivative of the convex function $\varphi$. The latter ability is especially useful when an algorithm stalls during the search for an optimum, and we can assess the proximity of the latest design from the algorithm to the optimum. However, while there are algorithms mathematically guaranteed to find the optimum, they tend to work for relatively simple models, and certainly do not work well for PK/PD models with a large number of covariates or when the objective function is non-differentiable, such as when we search for a minimax or standardized maximin design [25,26]. Details, examples and some technical justifications for the elegant theory are available in optimal design monographs, such as [27,28,29], among many others.

In contrast, no unified theory exists for finding exact optimal designs. In the rare instances where it is possible to derive an optimal exact design, the proof always depends sensitively on every aspect of the model assumptions, and the proof of optimality of a design for one model does not extend to another model, even if the model is only slightly changed. As far as we know, there are no algorithms that can guarantee they will find the optimal exact design. This is especially true for complicated nonlinear models for PK studies, where invariably, there are random effects in the model, so they have to be integrated out from the information matrix, implying that there are no closed forms for the matrix. The information matrix has to be approximated accurately using various ways, and computation issues become complex; these are well discussed in [10,16,17,18]. Further, the PK/PD studies are longitudinal with corrected errors, and there are usually multiple restrictions imposed on the study design. For example, there may be a maximum number of sampling points allowed per subject, the subjects have to be grouped according to different sampling schedules for evaluation, and sometimes doses or sampling times cannot be too close to one another. The problem is no longer a convex optimization problem, and new effective tools need to be developed to find efficient designs and analyze data for such studies. As far as we know, there are no free statistical software packages that fulfill such tasks.

PopED is a most timely and free package with many capabilities, specially designed for solving design and analysis issues in PK/PD studies. It is a continually developing package with state-of-the-art techniques and widely used by pharmacometricians but not in the statistical and mathematical research communities. It uses the optimization algorithm that evaluates the design based on the Fisher information matrix of the PK model structure and the range of the parameters of the model. For the rest of this section, we describe the functionality of the optimization tool of the R library PopED to build and evaluate PK models. To this end, we provide a general guide on some feature commands in the PopED library that were used in the analyses.

Building a population PK model in PopED starts with model definition. Model definition involves four functions: the model-defining function, parameter-defining function, residual-defining function, and the PopED database-building function. The model-defining function takes the values of $t_{ij}$ as variables from the PopED database, and returns $Y_{ij}$ as the outcome. The parameter function specifies the PK between-individual parameter (${\theta }_i$) and the within-individual parameter ($U_i$) to be used for the simulations. This function defines ${\theta }_i$ as a mathematical function composed of the covariate $A_i$ and the random effects component ${\eta }_i$. The residual function adds residual variability ${ϵ}_{ij}$ to the model, and both additive and proportional residual variabilities are allowed.

The last and important component is the PopED database function. This function is an all-encompassing function, where the PopED database is created by specifying the values to be used for the simulation. The PopED database function specifies design-related settings that include the initial values for the mean parameters. The values for the mean between-individual parameter ($\overline{\theta }$), random effects variances ($\mathbf{G}$), and the residual variability (${\sigma }^2$) are set for the simulation, along with the initial number of arms (groups), group size, doses, dose intervals, and sample times as part of the design specifications. If the design optimization is to be performed as part of the simulation, the range of doses and the dose interval, and potential sampling window per sampling points are also set at this stage.

3.2. Evaluation and Optimization Metric

Based on the pre-designated mean or median PK parameters and other dose- and time-related specifications defined in the PopED database, the simulation function generates a concentration–time plot with 95% prediction confidence bands for the model. With this initial simulation performed, PopED provides a design evaluation feature for comparing various designs. The evaluation function returns three evaluation criteria for each design: its objective function value (OFV), its Fisher information matrix (FIM), and the relative standard error (RSE) for model parameters. In our case, OFV is the log value of the determinant of the FIM, and RSE is the standard error of the estimated parameters as a percentage of the parameter estimate, i.e., $RSE\%\left(x\right)=\left(SE\right(x)/x)·100$.

The RSE(%) values are reported for the estimated fixed-effect and random-effect parameters. Larger OFV values are associated with better designs, and smaller RSEs for parameters imply that estimates are more accurate, so the design is more efficient than designs with a larger value. When p is the number of PK parameters in the model, PopED has a function to compare the efficiency of an alternative design relative to the initial design using the measure

$$efficiency={\{\frac{exp\left({OFV}_{alt}\right)}{exp\left({OFV}_{init}\right)}\}}^{\frac{1}{p}}.$$

(4)

In practice, efficiency denotes the relative size of the sample needed for the alternative design against the initial design, indicating the number of times the alternative design has to be replicated to achieve the same performance level as the initial design.

PopED offers a determinant-based method for continuous and discrete optimization, i.e., it assumes that the design criterion is of D-optimality and the sought design maximizes the determinant of the information matrix among all designs. The design optimization function provides the sampling times that maximize the OFV of the initial design. This means that the time points provided by the optimization function are those that maximize the information content of the trial. Selecting those points as sampling times will maximize the precision of the estimates for the model parameters, as the OFV value is proportional to the determinant of FIM. Depending on the options, the optimization can be over the choice of dose levels or dose intervals. Similarly, these options provide the best dose regimen for the maximized information content and hence the most precise statistical inference for the model parameters.

4. Results

We use the two applications from the previous section to illustrate how we use PopED to compute the optimal designs and evaluate the relative efficiencies of other designs.

4.1. Application 1: One-Compartment First-Order Absorption Single and Multiple Dosing PK Model

Consider a one-compartment model with first-order absorption and no covariate. The nominal values we used in the simulation for the fixed-effect parameter means, the random effects variances, and the residual variability variance are borrowed from the introductory website of PopED and shown below:

Fixed effects parameter mean: ($\overline{V},{\overline{K}}_A,\overline{CL},\overline{F}$)	(72.8, 0.25, 3.75, 0.9)
Random effects variance (G)	$diag(0.09,0.09,0.0625)$
Residual variability variance (${\sigma }^2$)	0.04

Suppose we have 40 patients, and we follow them five times for 250 h after the initial dose of 20 mmol on the dose interval of 24 h. The design problem is to select five optimal sampling time points from the time range $[0,250]$ in hours after drug administration. We consider two designs: one with well-dispersed sampling points over the observation period, and the other with a skewed distribution of the sampling time points. A well-dispersed design has spread out sampling time points after the initial dose. An example of such a design is one with five equally allocated points at 50, 100, 150, 200, and 250 h after drug administration. An example of a design with a skewed distribution of sampling time points is the design that takes observations at 6, 12, 24, 216, and 240 h after the initial dose. We evaluate the relative efficiencies of the two designs and examine the optimized sampling times with and without the time condition per each sample. Figure 1 displays the simulated plot of the dose–concentration of the initial design, where the shaded area denotes the $95\%$ confidence band.

Table 1. Sampling time optimization result (Application 1).

Design	Sampling Times	OFV	Efficiency
Design 1 (dispersed)	(50, 100, 150, 200, 250)	24.02	1
Design 2 (skewed)	(6, 12, 24, 216, 240)	37.49	0.15
Design 3 (optimized)	(0.57, 14.15, 14.15, 14.15, 240)	40.79	0.09
Design 4 (optimized with condition)	(15.31, 51.02, 146.9, 168, 242.9)	38.11	0.13

summarizes the simulation results of the four experimental designs for the model. Design 1 and Design 2 are respectively well dispersed and skewed as explained previously. For Design 1 and Design 2, no other condition was imposed on the sampling window for each point. All five sampling points were free to be selected from the range of the initial design space [0, 250]. The OFV value for the dispersed design (Design 1) is 24.02, whereas for the more skewed design (Design 2), it is 37.49. This suggests that the skewed design provides more information for the same number of sampling points and participants. The efficiency is 0.15, which suggests that for the alternative skewed design to provide the same informativeness as the initial (dispersed) design, the alternative design would only have to have 15% of the participants in the initial design for the trial. The result highly suggests that for this experiment, the design with sampling points concentrated toward the beginning and the end of the experiment performs better.

We can similarly calculate the optimized sampling points for this design. Design 3 shows the optimized sampling points for this model, which are 0.57, 14.15, 14.15, 14.15, and 240 h. These sampling points yield the OFV of 40.79 with an efficiency of 0.09, meaning that theoretically, the experiment would only need 9% of the initial study participants if the experiment uses the optimized sampling points. However, having duplicate sampling time points—in the case of Design 3, the three duplicates of 14.15—is problematic. To resolve this matter, we can think of a design optimization with some restrictions on the sampling windows for each point. For Design 4, the condition was added that the five sampling points should come from an equally divided and minimally overlapping subspace of the initial design space. The first to fifth sampling points were bound to be selected from [0, 50], [50, 100], [100, 150], [150, 200], and [200, 250], respectively. The optimized sampling points with such a condition were 15.31, 51.02, 146.9, 168, and 242.9 h with OFV of 38.11 and efficiency of 0.13 (Design 4). Although the OFV was lower in the case of optimization with the condition than the unconditioned optimization, Design 4 provides more intuitive and feasible sampling points with a relatively small difference in efficiency.

Once the optimal sampling time points are identified, further analysis of the optimal dose regimen can be performed.

Table 2. Optimized dose regimen for (Application 1).

Design	Sampling Time	Dose/Dose Interval	OFV	Efficiency
Design 4	(15.31, 51.02, 146.9, 168, 242.9)	20 mmog/24 h	38.11	1
Design 4–1		40 mmog/48 h	43.45	0.47
Design 4–2		30 mmog/36 h	42.50	0.53
Design 4–3		10 mmog/12 h	36.12	1.33
Design 4–4		5 mmog/6 h	33.35	1.97

shows the design evaluation results for different dose regimens. Designs 4–1 to Design 4–4 have the same sampling points as Design 4 and differ only in the doses and dose intervals. Different variations of dose and dose interval combinations reveal that in our model, when smaller doses of the drug are administered more frequently, the design performs increasingly worse, and as the dose administration becomes less frequent with a larger dose, its performance improves. Design 4–1 is the best dose regimen out of the four variations in terms of efficiency, as it requires only 47% of the sample size as Design 4 to provide the same information.

4.2. Application 2: Single Dosing Two-Compartment PK Model with Perturbation

Shotwell (2016) used a two-compartment PK model with perturbation and a single dosing scheme to study the effect of piperacillin on hemodialysis by injecting the antibiotic into patients with reduced kidney function [24]. For this simulation, we assume an eight-hour dosing cycle with 100 participants. The cycle starts with an initial administration of the drug by intravenous infusion that lasts for 60 min at a dose of 3000 mg/h to the central compartment. A 120 min dialytic clearance (perturbation) takes place from 180 to 300 min (3 to 5 h) from the initial administration at 15 L/h rate ($k_D$). The nominal values used in the simulation assume that the target population is patients receiving hemodialysis, and are from Shotwell (2016) shown below [24]:

Fixed effects parameter median: (${\tilde{v}}_1,{\tilde{k}}_{10}{\tilde{k}}_{12},{\tilde{k}}_{21}$)	(10, 0.35, 3.5, 1.5)
Random effects variance (G)	$diag$(0.06, 0.06, 0.06, 0.06)
Residual variability variance (${\sigma }^2$)	0.04

Shotwell (2016) suggested that the vector of the fixed-effect parameter median for patients with renal insufficiency be 10, 0.35, 3.5, and 1.5. In their research, the coefficient of variation for the fixed-effect parameter was set to 25%. The relationship between the coefficient of variation ($CV(\%)$) and the variance of random effects for inter-individual variability ($var\left({\eta }_i\right)$) is $CV(\%)=\sqrt{exp\left(var\left({\eta }_i\right)\right)-1}·100(\%)$ [30]. Using this relationship, the random effects variances were all set to 0.06, and the residual variability variance was again set to 0.04.

The simulation picks 2, 4, 6, and 8 h after the initial dose as the sampling times (initial design). The sampling times for the alternative design are at 0.5, 1, 6, and 8 h, implying more samples near the start of the dosing cycle and around hemodialysis. Figure 2 and Figure 3 display the concentration plot of the central and peripheral compartment for the initial design, respectively, and the shaded area in each of the plots corresponds to a 95% confidence band.

Table 3. Sampling time optimized results (Application 2).

Compartment	Design	Sampling Time	OFV	Efficiency
central	initial	(2, 4, 6, 8)	49.12	1
	alternative	(0.5, 1, 6, 8)	47.31	1.22
	optimized (continuous)	(0.0008, 1.13, 3.46, 8)	60.38	0.29
	optimized (discrete)	(1, 2, 4, 8)	53.70	0.60
peripheral	initial	(2, 4, 6, 8)	46.93	1
	alternative	(0.5, 1, 6, 8)	52.72	0.53
	optimized (continuous)	(0.01, 3.04, 5.13, 8)	58.30	0.39
	optimized (discrete)	(1, 3, 5, 8)	54.09	0.55

lists the optimized sampling times, along with the OFV, and efficiency values for the various designs. For both compartments, the initial design has a set of well-dispersed sampling times at 2, 4, 6, and 8 h. In contrast, the alternative design has more concentrated sampling time points. The alternative design for the central compartment has an efficiency of 1.22, implying that the initial design outperforms the alternative design by 22%. For the peripheral compartment, the opposite result is true; for it to have the same performance as the initial design, the alternative design requires only 53% of the number of participants in the initial design.

The continuous optimized sampling points for the central compartment are at 0.0008, 1.13, 3.46, and 8 h, and the corresponding sampling time points using discrete optimization are at 1, 2, 4, and 8 h. For both optimization schemes, the resulting efficiencies are 29% and 60%, respectively, showing improvement. For the peripheral compartment, the optimized sampling times are at 0.01, 3.04, 5.13, and 8 h for the continuous optimization, and are at 1, 3, 5, and 8 h using discrete optimization. This shows that optimization for the peripheral compartment has resulted in a 39% and a 55% efficiency improvement relative to the initial design.

For the central compartment, the continuous and discrete optimized sampling points clearly show enhanced performance compared to both the initial and alternative designs. However, we observe that in the simulation for the peripheral compartment, the differences in performance between the alternative design and the optimized designs are more nuanced, making it a viable option to choose the alternative design if the sampling times suggested by the two optimizations, while being optimized in terms of performance, is simply infeasible or impractical.

Table 4. Relative standard error (RSE%) of the estimated parameters (Application 2).

Compartment	Design	${\mathit{k}}_{10}$	${\mathit{k}}_{12}$	${\mathit{k}}_{21}$	${\mathit{v}}_1$	$\mathit{var}\left({\mathit{\eta }}_1\right)$	$\mathit{var}\left({\mathit{\eta }}_2\right)$	$\mathit{var}\left({\mathit{\eta }}_3\right)$	$\mathit{var}\left({\mathit{\eta }}_4\right)$
central	initial	164	257	38	161	37	133	77	31
	alternative	68	227	136	94	55	127	248	29
	optimized (continuous)	4	7	7	3	35	57	56	20
	optimized (discrete)	55	97	21	54	37	74	53	26
peripheral	initial	24	46	37	5	311	240	169	25
	alternative	23	16	14	7	277	54	71	24
	optimized (continuous)	4	23	19	3	32	209	139	20
	optimized (discrete)	84	149	32	82	35	77	66	26

displays the RSE(%) per fixed- and random-effect parameter. It shows the proportion of the standard error of each parameter compared to its estimate in percentage. Generally, the RSE for continuous optimized design outperforms discrete designs. This is to be expected, as the continuous design space allows for the more flexible selection of the sampling and dosing regimen. Discrete optimization is performed in limited design space, which leads to it being less precise than when the continuous design space is used. Overall, RSE is dropped significantly after optimization for both compartments, except for peripheral discrete optimization. Again, the drop in RSE is more significant in the central compartment than it is in the peripheral compartment. This suggests that if the purpose of the study is to investigate the PK relationship in the peripheral compartment, choosing the alternative design could serve this purpose if the sampling windows do not allow for the optimized sampling times.

5. Discussion and Conclusions

This article demonstrated the utility of non-linear mixed effects models for studying population pharmacokinetics with two applications. The first is illustrative using a one-compartment first-order PK model with a single oral dosing and multiple bolus administration. The second is a new application for a single dosing two-compartment PK model with perturbation. The simulation identified the best design with the optimal sampling points or dose regimen based on the objective function value, along with efficiencies of the initial and optimized designs. The use of the design evaluation and optimization tool PopED is presented in the R setting, the statistical programming language that is now widely used even in clinical trials and experimental designs. An efficient computing tool such as PopED is integral in drug development in finding which design is the most efficient given the available resources in the trial.

Overall, the evaluation and optimization for various population PK models are facilitated by the PopED library in R in that it provides a platform for the users to specify flexible models for a PK study, and has tools to design and analyze the data. The library operates in an open-source environment R, and this further adds flexibility to work with other statistical libraries in R for a more in-depth analysis. However, this also means that the familiarity of the programming language R is a requirement to fully use PopED, and this may limit its user base to R users. In contrast, the NLMEM Population PK/PD modeling software provides its own user interface, which eliminates the need to have proficiency in a certain programming language.

We note that the library does not provide estimates of the parameters for the population PK/PD model from a dataset. If parameter estimates from previous studies are not available, then an additional tool for model fitting and parameter estimation is needed. However, with the recent surge in demand for a good study optimization tool for cost-efficient study designs, software that specializes in optimization such as PopED is becoming more widely used in that it provides various optimization and simulation sub-analysis features, such as the analysis for irregular dosing, below the limit of quantification (BLQ) data analysis, power calculation, optimization of subjects per group, shrinkage, and Bayesian FIM estimates. In this sense, the population PK design evaluation and optimization using PopED is both a helpful and useful tool for performing population PK design analysis.

Future work includes extending the current features in PopED to perform additional tasks, such as power calculation or finding Bayesian or minimax types of optimal designs. Further enhancements can include a package that performs a comprehensive sensitivity analysis to ascertain the design sensitivities to all aspects of the model assumptions.

We close by noting that while the optimal sampling design is a valuable tool to assess the study design based on FIM, the current gold standard approach is to use stochastic simulation and estimation (SSE) to assess the PK/PD study design. The SSE allows researchers to assess both the bias and precision of parameters respective to the study design, where the optimal sampling design only allows precision assessment. Future work should also provide insights into the two approaches and how to integrate them more effectively.

For now, we believe that both approaches are valid and useful, and they should be used jointly for research. First, the optimal sampling design is derived mathematically or using an algorithm based on a criterion that may not fully incorporate the goals of the study in practice. They are invariably model based, and so optimal designs that are theoretically derived based on mathematics alone may not be realistic to implement. For example, they may have too few sampling time points than the number that clinicians like. Optimal designs therefore should be used only as a guide or as a rough benchmark to assess other designs. In practice, the implemented design should meet practical needs and stray not too far from the optimal design; otherwise, the implemented design can lose substantial statistical efficiency. When there are several implementable designs proposed by clinicians, or otherwise, one may select the design with the highest efficiency relative to the optimum so that the precision of the estimates is the highest to the extent possible.