In john-harrold/ubiquity: PKPD, PBPK, and Systems Pharmacology Modeling Tools

knitr::opts_chunk$set(echo = TRUE, message=FALSE, eval=FALSE)
require(ubiquity)
#require(deSolve)
require(ggplot2)
#require(foreach)
#require(doParallel)
#require(rhandsontable)

Introduction

Sometimes we need to analyze data where the independent variable is not time. You can do this with ubiquity and this tutorial will highlight how that is done. If you have not done so already, please review the naive-pooled parameter estimation tutorial because this will build on the concepts covered there. The workshop (workshop.ubiquity.tools) provides an example of how to analyize static in vitro data. To make a copy of these scripts and other supporting files in the current working directory run the following:

library(ubiquity)
fr = workshop_fetch(section="In Vitro", overwrite=TRUE)

This should create the following files in your current working direcotry:

• system-in_vitro.txt - System file that can be run both in vivo and in vitro
• in_vitro_er_data.csv- In vitro data file
• mk_data_in_vitro.R- Script to generate in vitro data
• analysis_in_vitro.R- Script to perform parameter estimation using the in vitro data

This example uses a system from Chakraborty and Jusko. J Pharm Sci 91(5), 1334-1342. Specifically equation 5 from that article where the effect of two drugs, A and B, is driven by their concentrations ($C_A$, and $C_B$):

$$ E = 100 \left( 1- \frac{ \frac{I_{max,A}(C_A) ^{\gamma_A}}{(\Psi IC_{50,A})^{\gamma_A}} + \frac{I_{max,B}(\xi C_B)^{\gamma_B}}{(\Psi IC_{50,B})^{\gamma_B}} \ + (I_{max,A} + I_{max,B} + I_{max,A} I_{max,B}) \ \times \frac{(C_A) ^{\gamma_A}}{(\Psi IC_{50,A})^{\gamma_A}} \times \frac{(\xi C_B)^{\gamma_B}}{(\Psi IC_{50,B})^{\gamma_B}} \ }{ \frac{(C_A) ^{\gamma_A}}{(\Psi IC_{50,A})^{\gamma_A}} + \frac{(\xi C_B)^{\gamma_B}}{(\Psi IC_{50,B})^{\gamma_B}} \ +\frac{(C_A) ^{\gamma_A}}{(\Psi IC_{50,A})^{\gamma_A}} \times \frac{(\xi C_B)^{\gamma_B}}{(\Psi IC_{50,B})^{\gamma_B}} +1 } \right)
$$

$$ \xi = \frac{IC_{50,A}}{IC_{50,B}} $$

This has been implemented in the PKPD system file (system-in_vitro.txt) shown at the bottom. This implementation is dynamic meaning the effect changes with the PK of the drug. However the effect is considered instantaneous and typically analyzed from in vitro data using algebraic relationships. For this system consider the following data:

er_data = read.csv(system.file("ubinc", "csv","in_vitro_er_data.csv", package="ubiquity"))
er_data$C_B0 = as.factor(er_data$C_B0)
p = ggplot(data=er_data) +
    geom_point(aes(x=C_A0, y=Effect, group=C_B0, color=C_B0)) +
    geom_line(aes(x=C_A0, y=ave_eff, group=C_B0, color=C_B0)) + 
    xlab("CA (ng/ml)") +
    scale_color_discrete(name="CB (ng/ml)") +
    theme_minimal()
p = gg_log10_xaxis(p)
p = prepare_figure(p, purpose="shiny")
print(p)

Here both the concentrations of drug A and B are altered (independent variable) and the Effect is measured (dependent variable).

The sysem file: system-in_vitro.txt

First lets discuss the way the system file is structured. There is bolus dosing for the PK specified, but the initial conditions of the effect compartments are also defined in terms of system parameters:

<I> Cp_A = C_A0
<I> Cp_B = C_B0

These system parameters (C_A0 and C_B0) have a default value of zero. So by default the system would run like any in vivo PKPD system. But these initial condition placeholders will be used when performing in vitro analyses. This allows you to use the same system file for both in vitro and in vivo analyses. This is useful because you do not have to do your static in vitro analysis in one system file and then copy and paste to your in vivo dynamic system.

The dataset: in_vitro_er_data.csv

This table contains a snapshot of the relevant columns of the dataset:

tmpdata =  head(er_data, n=10)
tmpdata[["ave_eff"]] = NULL

flextable::flextable(tmpdata) |> flextable::autofit()

The C_A0 and C_B0 columns correspond to the concentrations that elicit the Effect. The treat column is a unique name for the combination of the two drug concentrations. It has only letters numbers and underscores. This is intentional so I can use this as a cohort name when I'm constructing the estimation script below. The samp_time column is set to 1. This is arbitrary because the mapping in the cohort definitions requires a time column. You'll see why it is arbitrary below.

The analysis script: analysis_in_vitro.R

This analysis scripts has some aspects that are unique to the in vitro analysis being performed.

Output times

First is the output times. When you run simulations (which is what happens when performing a parameter estimation) more output times results in slower simulations. So it is important to only include the ncesssary output times. In this case we only have two output times (0 and 1):

cfg=system_set_option(cfg, group  = "simulation", 
                           option = "output_times", 
                           seq(0,1,1))

The final value of 1 was chosen because that corresponds to the time column (samp_time) in the analysis dataset.

Making a dynamic simulation static:

Next we set the simulation option dynamic to FALSE:

cfg=system_set_option(cfg, group  = "simulation", 
                           option = "dynamic", 
                           value  = FALSE)

This will fix the values of the differential equations to 0 for the purposes of the subsequent simulations. Only available in ubiquity version $\ge$ 2.04.

Datasets

For the dataset we are reading it into a dataframe and loading it that way. This was done so the dataset could be used internally with ubiqiuty and also to construct the analysis below:

er_data = readr::read_csv("in_vitro_er_data.csv")
cfg = system_load_data(cfg, dsname     = "er_data",
                            data_file  = er_data)

Defining cohorts

In the estimation tutorial cohorts are defined individually. We could do that here too but it would be rather tedious. See we need a cohort for every unique combination of C_A0 and C_B0 in the dataset. To do this we are going to loop through each unique value in the treat column and create a cohort. The variable tmp_treat contains the value of the current treatment so the first thing we do is get the records for the current treatment:

treat_recs = dplyr::filter(er_data, treat == tmp_treat)

This subset of the data is used to define the initial condition parameters for the current cohort using the cp field (see the help for system_define_cohort() for more information about this option). Because of the way the treat column was constructed we can use this as the cohort name. This will allow us to link the simulated output below to the original dataset in the post processing section below. The inputs are set to NULL here and we do not have to change them because there is no dosing.

cohort = list(
    name         = tmp_treat,
    cf           = list(
                     treat     = c(tmp_treat)),
    cp           = list(
                     C_A0 = treat_recs$C_A0[1],
                     C_B0 = treat_recs$C_B0[1]),
    inputs       = NULL,
    outputs      = NULL,
    dataset      = "er_data")

Because we set the simulation option dynamic to FALSE, the initial condition is set to the values for the current treatment, and there are no inputs, then the system is effectively behaving like an in vitro system.

Postprocessing

You cannot rely on the normal figure generation and reporting elements. In this example we can take the simulated output at the estimation erp and the original data dataset er_data to create meaningful diagnostic plots. This is because we can link the two datasets using the treat column from the original dataset to the COHORT column in erp$pred:

df_orig = er_data |>
  dplyr::select(C_A0, C_B0, ave_eff, treat) |>
  dplyr::distinct()

df_est = erp$pred |>
  dplyr::filter(!SMOOTH) |>
  dplyr::rename(treat = COHORT)

df_plot = dplyr::full_join(df_est, df_orig, by="treat") |>
  dplyr::mutate(C_B0 = as.factor(C_B0))

library(ggplot2)
p = ggplot(data=df_plot) +
    geom_point(aes(x=C_A0, y=OBS,  group=C_B0, color=C_B0)) +
     geom_line(aes(x=C_A0, y=PRED, group=C_B0, color=C_B0)) + 
     scale_x_log10()

Contents of system-in_vitro.txt

cat(readLines(system.file("ubinc", "systems","system-in_vitro.txt", package="ubiquity")), sep="\n")

john-harrold/ubiquity documentation built on March 1, 2025, 12:35 p.m.