PKPlus

library(tidyverse)
library(gastroPlusAPI)

if (rstudioapi::isAvailable()){
  current_working_directory <- dirname(rstudioapi::getSourceEditorContext()$path)
  setwd(current_working_directory)
}

gpx_service <- start_service(verbose=FALSE)

✔ Configured the GastroPlus Service

gpx_service$is_alive()

[1] TRUE

project_path = "../../ProjectFiles/Midazolam.gpproject"

input_names <- c("5mg_iv", "30mg_sol", "7.5mg_sol")
observed_data_group_names <- c("5 mg iv", "30 mg Solution", "7.5 mg Solution")
observed_data_series_names <- c("Plasma Conc", "Plasma Conc", "Plasma Conc")
dose_type <- c(PKPlusDoseRoute$Intravenous, PKPlusDoseRoute$Oral, PKPlusDoseRoute$Oral)
dose_mg <- c(5, 30, 7.5)
infusion_time_h <- c(0, 0, 0)
body_weight_kg <- c(70, 70, 70)

pkplus_run_name <- "Run_for_1IV_2PO"
use_individual_BW <- TRUE
fit_F_for_oral_doses <- TRUE
fit_lag_time_for_oral_doses <- TRUE
fraction_unbound_plasma <- 0
elimination_model_type <- PKPlusEliminationModel$Linear
analyses_methods <- c(PKPlusAnalysisMethod$NonCompartmentalAnalysis,
                      PKPlusAnalysisMethod$OneCompartment,
                      PKPlusAnalysisMethod$TwoCompartment,
                      PKPlusAnalysisMethod$ThreeCompartment)

#Load the project
open_project(project_path)

#Loop begins for PKPlus inputs
for (i in (1:length(input_names))){

  #Find the assigned observed data
  tryCatch({
    exposure_series <- get_experimental_data_inventory(group_type = ObservedDataGroupType$ExposureData) %>%
      filter(
        group_name == observed_data_group_names[i],
        series_name == observed_data_series_names[i]
      )
  }, error = function(e) {
    message(
      "The assigned observed data is not found, please check if it truly exists or possible typos."
    )
  })
  
  #exposure_series should contain a single row, check if it truly is
  if (nrow(exposure_series) > 1) {
    cli::cli_abort("exposure_series contains more than one data, please check!")
  } 
  
  #Assign the PKPlus input details
  dose_amount <- ScalarObservable$new(
    scalar_type = ScalarType$Mass,
    scalar_value = ScalarValue$new(value = dose_mg[i], unit = "milligram")
  )
  body_weight <- ScalarObservable$new(
    scalar_type = ScalarType$Mass,
    scalar_value = ScalarValue$new(value = body_weight_kg[i], unit = "kilogram")
  )
  infusion_time <- ScalarObservable$new(
    scalar_type = ScalarType$Time,
    scalar_value = ScalarValue$new(value = infusion_time_h[i], unit = "hour")
  )
  
  series_key <- SeriesKey$new(
    group_type = exposure_series$group_type,
    group_name = exposure_series$group_name,
    series_type = exposure_series$series_type,
    series_name = exposure_series$series_name
  )
  
  
  input_config <- PKPlusInputConfiguration$new(
    dose_route = dose_type[i],
    exposure_series_key = series_key,
    dose_amount = dose_amount,
    infusion_time = infusion_time,
    body_weight = body_weight
  )
  
  add_pkplus_input(input_names[i], input_config)
  
  #Info about the loop to user
  cli::cli_alert_info(
    "Input added: {input_names[[i]]}, Dose Amount: {dose_mg[i]}, Body Weight: {body_weight_kg[i]}, Infusion Time: {infusion_time_h[i]}"
  )
  
} #Loop ends

ℹ Input added: 5mg_iv, Dose Amount: 5, Body Weight: 70, Infusion Time: 0

ℹ Input added: 30mg_sol, Dose Amount: 30, Body Weight: 70, Infusion Time: 0

ℹ Input added: 7.5mg_sol, Dose Amount: 7.5, Body Weight: 70, Infusion Time: 0

#Check if inputs implemented successfully
if( length(get_pkplus_inputs()$pkplus_input_name) != length(input_names)){
  cli::cli_abort("An error occurred during input, please check your inputs.")
}


#Create PKPlus Run and Execute
add_pkplus_run(pkplus_run_name = pkplus_run_name, pkplus_input_names = input_names)

run_config <- PKPlusRunConfiguration$new(use_individual_body_weights = use_individual_BW,
                                          fit_bioavailability = fit_F_for_oral_doses,
                                          fit_time_lag = fit_lag_time_for_oral_doses,
                                          fraction_unbound_plasma = fraction_unbound_plasma,
                                          kinetics_type = elimination_model_type,
                                          model_types = analyses_methods)

configure_pkplus_run(pkplus_run_name = pkplus_run_name, pkplus_run_config = run_config)

tryCatch({
  execute_pkplus_run(pkplus_run_name = pkplus_run_name)
  cli::cli_alert_info("Run executed: {pkplus_run_name}, Inputs: {input_names}")
}, error = function(e) {
  message("A problem occurred while executing PKPlus run: {pkplus_run_name}")
})

ℹ Run executed: Run_for_1IV_2PO, Inputs: 5mg_iv, 30mg_sol, and 7.5mg_sol

pkplus_output <- get_pkplus_run_output(pkplus_run_name = pkplus_run_name)

non_compartmental_analysis_output <- pkplus_output$non_compartmental_analysis_output

#Loop_1 for data manipulation of NCA outpus
for(i in (1:length(input_names))){
  
NCA_data_longer <- as.data.frame(non_compartmental_analysis_output[[i]]) %>%
  mutate(across(everything(), as.character)) %>%
  pivot_longer(
    cols = everything(),
    names_to  = "column1", 
    values_to = "column2"
  ) %>% filter(column1 != "dose_route") #delete dose_route, unnecessary and violating data manipulation

#Warning occurrs for the following code due to NA assignment, masking it.
NCA_df <- suppressWarnings(
  NCA_data_longer %>%
    separate(
      column1, 
      into = c("parameter", "attribute"), 
      sep = "\\.(?=[^.]+$)"   # split at the last dot
    ) %>%
    pivot_wider(
      names_from  = attribute,
      values_from = column2
    ) %>%
    select(parameter, value, unit)
)

NCA_df_name <- paste0("NCA_", i)
assign(NCA_df_name, NCA_df)

cli::cli_alert_info("NCA output is generated: {NCA_df_name}, Inputs: {NCA_df[1,1]}")
} #Loop_1 ends

ℹ NCA output is generated: NCA_1, Inputs: X.30mg_sol.

ℹ NCA output is generated: NCA_2, Inputs: X.5mg_iv.

ℹ NCA output is generated: NCA_3, Inputs: X.7.5mg_sol.

#Loop_2 for combining NCA outputs and copying to clipboard
NCA_df_list <- list()
for (i in (1:length(input_names))) {
NCA_df_list[[i]] <- get(paste0("NCA_", i))
} #Loop_2 ends

combined_NCA_df <- do.call(rbind, NCA_df_list)

write.table(
  combined_NCA_df, 
  file      = "clipboard", 
  sep       = "\t", 
  row.names = FALSE)
cli::cli_alert_info("NCA outputs were copied to clipboard")

ℹ NCA outputs were copied to clipboard

glimpse(combined_NCA_df)

Rows: 42
Columns: 3
$ parameter <chr> "X.30mg_sol.", "auc_first_moment_last", "auc_infinity", "auc…
$ value     <chr> NA, "0.00221953", "0.000429704", "0.00040325", "70", "69815.…
$ unit      <chr> NA, "(mg/mL)*h*h", "(mg/mL)*h", "(mg/mL)*h", "kg", "mL/h", "…

observed_conc_series <- pkplus_output$compartmental_analysis_output$observed_concentration_series

one_CA_series <- pkplus_output[["compartmental_analysis_output"]][["model_outputs"]][[1]][[2]][["predicted_concentration_series"]]
two_CA_series <- pkplus_output[["compartmental_analysis_output"]][["model_outputs"]][[2]][[2]][["predicted_concentration_series"]]
three_CA_series <- pkplus_output[["compartmental_analysis_output"]][["model_outputs"]][[3]][[2]][["predicted_concentration_series"]]


p <- ggplot()

#Plotting observed data
for (i in (1:(nrow(observed_conc_series)))) {
  
  df_data <- observed_conc_series$series[[i]]
  
  df_data$Legend <- observed_conc_series$input[i]
  
  p <- p + geom_point(
    data = df_data,
    mapping = aes(
      x     = independent,
      y     = dependent,
      color = Legend
    ),
    size = 3.5, 
    alpha=0.75, 
    shape=22,
  )
}

#Plotting the best CA
if(pkplus_output[["best_model"]] == "OneCompartment") {
  best_CA <- "one_CA_series"
} else if(pkplus_output[["best_model"]] == "TwoCompartment") {
  best_CA <- "two_CA_series"
} else if(pkplus_output[["best_model"]] == "ThreeCompartment") {
  best_CA <- "three_CA_series"
}

for (i in (1:(nrow(get(best_CA))))) {
  
  df_data_CA <- get(best_CA)$series[[i]]
  
  df_data_CA$Legend <- get(best_CA)$input[i]
  
  p <- p + geom_line(
    data = df_data_CA,
    mapping = aes(
      x     = independent,
      y     = dependent,
      color = Legend
    ),
    linewidth=1.25
  )
}

#Styling
p <- p + theme_classic() + 
  labs(title= paste0("Best Model is ", pkplus_output[["best_model"]]), x= paste("Time"), y=paste("Concentration")) +
  theme(
    axis.title.x = element_text(face = "bold", size = 12),
    axis.title.y = element_text(face = "bold", size = 12),
    axis.text.x = element_text(face = "bold", size = 11),
    axis.text.y = element_text(face = "bold", size = 11),
    plot.title = element_text(hjust = 0.5, face = "bold", size = 12, color="#051082")
  ) 

ylog_p <- p + scale_y_log10()

p

ylog_p

Warning in scale_y_log10(): log-10 transformation introduced infinite values.
log-10 transformation introduced infinite values.
log-10 transformation introduced infinite values.
log-10 transformation introduced infinite values.

gpx_service$kill()

[1] TRUE