Overview

PKPDindex() provides a tool for fitting Emax models to PKPD data. This package examines the relationship between drug concentration and biological response by applying a set of predefined Emax models to the provided dataset. It allows customization in terms of model selection, parameter initialization, and plot settings.

Software requirements

PKPDindex() is implemented in R and available through CRAN. Users require an R installation and PKPDindex package to perform the analyses described in this protocol. The package can be installed directly from https://CRAN.R-project.org/package=PKPDindex.

Input requirements and data structure

The PKPDindex() requires a dataset containing PKPD indices and a corresponding response variable that reflects the antimicrobial effect. Data should be provided as a standard R data frame, with each row representing an experimental condition or individual observation.

The PKPDindex() converts all column names to lowercase to prevent case-sensitivity issues. It also checks the presence of both the y_column (default name: “response”) and the predictor columns x_columns (default column name: auc_mic, cmax_mic, t_mic). If the dataset uses different column names, users can specify them manually within the function.

The y_column should be the log10-transformed change in cfu/mL (Δlog10 cfu/mL). Alternatively, users may supply raw cfu/mL counts at the initial CFU_init and 24-hour CFU_24 time points, in which case the function will automatically calculate the Δlog10 cfu/mL.

The user must calculate the x_columns from their PK model before analysis, since the function does not generate them. These include area under the concentration–time curve divided by the MIC auc_mic, peak drug concentration divided by the MIC cmax_mic, and time the drug concentration stays above the MIC t_mic.

The package validates that the required columns are present and contain numeric data. It does not automatically screen for zero values, skewness, or outliers. Model fitting failures are captured and reported in the output.

Emax model variants

The PKPDindex() package defines eight nonlinear Emax models (m1–m8) (Table 1), each representing different combinations of key parameters: E0, Emax, EI50, and Hill coefficient. These models vary based on whether E0 and Emax are fixed or estimated during the fitting process. Additionally, some models incorporate the Hill coefficient to account for sigmoidal (S-shaped) dose-response curves. This flexibility allows the selection of the most appropriate model based on the complexity of the data and the biological system under study.

The eight Emax model variants represent alternative assumptions regarding baseline effect (E0), maximum effect (Emax), and the shape of the exposure–response relationship (Hill coefficient). In practice, there is no universally preferred model, and model selection should be driven by the observed data rather than predefined biological assumptions. Simpler models may be adequate when data are limited or do not support estimation of additional parameters, whereas more complex models may better capture nonlinear exposure-response relationships when sufficient information is available. The PKPDindex() package therefore fits all model variants and uses objective model-selection criteria to identify the most appropriate model for a given dataset.

Parameter initialization

Parameter initialization sets the starting values for Emax model fitting. By default, the Hill coefficient (γ) is initialized to 1, corresponding to the standard hyperbolic Emax model. This provides a neutral starting value for model fitting. E0 and Emax depend on the characteristics of the specific dataset. Therefore, users are required to specify initial values for E0 and Emax (E0_fix, Emax_fix) based on prior knowledge and data characteristics. Once these values are provided, the package automatically estimates an approximate EI50 starting value from the observed data. Alternatively, users may supply initial EI50 values through the EI50_init argument.

The EI50 estimate is derived from the PK/PD index value associated with a response halfway between the baseline effect (E0) and the maximum effect (Emax), calculated as:

Step 1: Prepare the dataset

Prepare a dataset containing the response variable and one or more PK/PD indices. By default, the response variable should be stored in a column named response and PK/PD indices in columns named auc_mic, cmax_mic, and t_mic. Alternative column names can be specified through the y_column and x_columns arguments. For demonstration purposes, the example dataset included with the package can be loaded using:

If different column names are used, they can be specified as follows:

Step 2: Specify initial parameter values

Before fitting the models, provide values for E0_fix and Emax_fix. Optional starting values for EI50 may also be supplied through the EI50_init argument. If EI50_init is omitted, the package will estimate starting values automatically. For example, the initial parameter values can be defined as follows:

Step 3: Fit Emax models

Execute the PKPDindex() function using the prepared dataset and specified initial parameter values. The package applies the nls() nonlinear least squares method in R to fit the models to the dataset. It iterates through each model, adjusts the parameters, and calculates the AIC and R² to assess the model fit.

The control settings for the NLS algorithm were configured as follows: The maximum number of iterations (maxiter = 500), the fitting tolerance (tol = 1e–5) and the minimum step size factor (minFactor = 1e–7). Users can modify these control settings to accommodate different data characteristics.

The package incorporates error handling to capture and report fitting failures, ensuring the analysis proceeds without interruption. Model fitting failures may occur when the nonlinear least squares algorithm is unable to identify parameter estimates that adequately describe the data, for example, because of insufficient information in the dataset, poor convergence, highly variable observations, or unsuitable starting values. Models that fail to converge are retained in the output and reported with a “Failure” status, while AIC, R², and parameter estimates are reported as NA. Such models are excluded from best-model selection. The current version of the package relies on the convergence assessment provided by the underlying nls() algorithm and does not perform formal sensitivity analyses of starting parameter values, apply parameter bounds or constraints, or conduct additional evaluations of local minima, parameter identifiability, or fit stability. Users may specify alternative starting values through the function arguments and should assess the biological plausibility of fitted parameter estimates when interpreting results (Algorithm 1).

Display full size

Algorithm 1

Example usage of the PKPDindex().

Step 4: Evaluate model performance

Model evaluation in the PKPDindex() package is based on two key performance metrics: the AIC and R². AIC is used for model comparison, as it balances goodness-of-fit with model complexity. A lower AIC value indicates a better model fit. R² is also calculated to provide an additional measure of how well the model explains the variability in the response data, with higher values indicating better explanatory power.

The model with the lowest AIC is selected as the best model for each PK/PD index. However, all model fits are summarized in the output, including their AIC, R², and success status, enabling transparent evaluation. The package also offers a model manual selection using the select_mod argument (Algorithm 2).

Display full size

Algorithm 2

Example usage of the PKPDindex().

Step 5: Generate plots

The function allows users to generate and customize the appearance of plots, including the colors of the data points point_color and fitted lines line_color, as well as adjust the axis limits (xlim and ylim), axis labels (x_label and y_label), and plot titles plot_title. Font sizes for titles title_cex, axis labels label_cex, axis tick labels axis_cex, and model annotation details detail_cex can also be adjusted. In addition, a log10 transformation can be applied to the x-axis for each PK/PD index through the argument log_scale_x to better visualize wide concentration ranges (Algorithm 3).

Display full size

Algorithm 3

Example usage of the PKPDindex().

Output structure

The PKPDindex() generates a list containing three key components. The All_Model_Results element provides a detailed summary of all models fitted for each PK/PD index, including performance metrics such as AIC, R², and parameter estimates, along with the fit status for each model. The Best_Models component summarizes the best-fitting model for each index, including the key parameters (E0, Emax, Hill coefficient, EI50) and corresponding evaluation metrics AIC and R² values. Lastly, the Plots element contains a list of plots for the best-fitting models. Table 2 shows an example of All_Model_Results output.

Visualization of results

When plot_results=TRUE, the PKPDindex() automatically generates plots that visualize the relationship between each PK/PD index (AUC/MIC, Cmax/MIC, T>MIC) and the pharmacodynamic response. Figure 1 shows an example of the expected plot output. For each index, the function displays the best-fitting model curve alongside the observed data points, allowing direct visual assessment of model performance. Each plot also displays the fitted model's key parameters and performance metrics including AIC, R², E0, Emax, γ, and EI50, which are annotated directly on the graph.

Display full size

Figure 1. Example of output plot from the PKPDindex(). The plot shows the optimal Emax models selected for each PK/PD index based on the lowest AIC values. Model performance is summarised by the corresponding R², E0, Emax, Hill coefficient, and EI50 estimates displayed on each plot.