Analysis of Static Systems, an In Vitro Example

Introduction

This tutorial requires ubiquity 2.04 or greater. Currently this is only available using the development version off of github

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:

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:

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.

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

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

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

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

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")

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

# Author: John Harrold <[email protected]>
#
# Analysis of an in vitro system. The system below was taken from:
#
# Chakraborty A, Jusko WJ. Pharmacodynamic interaction of recombinant human
# interleukin-10 and prednisolone using in vitro whole blood lymphocyte
# proliferation. J Pharm Sci. 2002 May;91(5):1334-42. doi: 10.1002/jps.3000.
# PMID: 11977109.
#

# #-------------#
# | Parameters  |
# #-------------#
#
# System parameters
#    name              value     lower      upper   units editable    grouping
#                                bound      bound
<P>  IC50_A              1.0       eps        Inf   -----      yes     Efficacy
<P>  IC50_B             50.0       eps        Inf   -----      yes     Efficacy
<P>  Imax_A              1.0       eps        Inf   -----      yes     Efficacy
<P>  Imax_B              1.0       eps        Inf   -----      yes     Efficacy
<P>  PSI                10.0       eps        Inf   -----      yes     Efficacy
<P>  G_A                 0.6       eps        Inf   -----      yes     Efficacy
<P>  G_B                 1.4       eps        Inf   -----      yes     Efficacy

<P>  Vp_A                1.0       eps        Inf   ml         yes     PK
<P>  CL_A                1.0       eps        Inf   ml/hr      yes     PK
<P>  Vp_B                1.0       eps        Inf   ml         yes     PK
<P>  CL_B                1.0       eps        Inf   ml/hr      yes     PK

<P>  C_A0                0.0       eps        Inf   -----      yes     IC
<P>  C_B0                0.0       eps        Inf   -----      yes     IC

# #-------------------#
# |Input Information |
# #-------------------#
#
# Bolus Events
# ------------
# times/events state   values    scale        units
<B:times>;              [  0  ];  1.0;        hours
<B:events>;      Cp_A;  [1.0  ];  1.0/Vp_A;   mg     
<B:events>;      Cp_B;  [1.0  ];  1.0/Vp_B;   mg     


# By default these values are zero but can be overwritten when 
# doing analysis of in vitro data.
<I> Cp_A = C_A0
<I> Cp_B = C_B0

<As> XI = IC50_A/IC50_B

<Ad> COMP_A    = SIMINT_POWER[Cp_A][G_A]/(SIMINT_POWER[PSI*IC50_A][G_A])
<Ad> COMP_B    = SIMINT_POWER[XI*Cp_B][G_B]/(SIMINT_POWER[PSI*IC50_A][G_B])
<Ad> COMP_Imax = Imax_A + Imax_B - Imax_A*Imax_B

<Ad> EFF_num    = Imax_A*COMP_A + Imax_B*COMP_B + COMP_Imax*COMP_A*COMP_B
<Ad> EFF_den    = COMP_A + COMP_B + COMP_A*COMP_B + 1.0
<Ad> EFF        = 100.0*(1.0-EFF_num/EFF_den)

# #-----------------------------#
# | ODEs, and State Information |
# #-----------------------------#

<ODE:Cp_A>  CL_A/Vp_A*Cp_A
<ODE:Cp_B>  CL_B/Vp_B*Cp_B

# #---------#
# | Outputs |
# #---------#
 
<O> Cp_A_mg_ml    = Cp_A
<O> Cp_B_mg_ml    = Cp_B
<O> Effect        = EFF

#<VP>  prop_err            1.0       eps        Inf   -----      yes     IC
#<OE:Effect> prop=prop_err

# #---------#
# | Options #
# #---------#
# specify different time scales
<TS:hours> 1.0