Plotting PK/PD hysteresis with variability in R using ggplot

The identification of hysteresis in a PK/PD relationship provides information on a possible delay between the plasma concentration and the effect.

The identification of hysteresis can further assist us in structural PK/PD model development by the inclusion of an effect compartment.

More information on hysteresis can be found on the following pages:

https://www.certara.com/2011/10/26/what-is-hysteresis-in-pkpd-analysis/?

https://www.ncbi.nlm.nih.gov/pubmed/24735761

In this post, I will show the R code required to generate some hysteresis figures, including the errorbars, using the ggplot2 package in R.

Dataset

I have created a completely random dataset of 5 individuals with the corresponding plasma concentrations (CP), a PD response which lags behind, and a direct effect. A screenshot of the dataset can be seen below and it can be downloaded for you to test.

R code

Lets load in the required libraries and explore the dataset with some standard figures in which we plot the concentration and the effect over time:

#Load in libraries needed in this script
library(ggplot2)
library(dplyr)
library(cowplot)

# Read the csv file
df <- read.csv('PKPD_dataset.csv')
df$ID <- as.character(df$ID) # ID needs to be a character for plotting

# Concentration over time with seperate lines per ID
p1 <- ggplot(df,aes(x=TIME,y=CP,color=ID))+
  geom_line()+
  xlab("Time after dose (h)")+
  ylab("Concentration")+
  theme_bw()+ 
  theme(legend.position = "none") 

# Effect over time (lag effect)
p2<- ggplot(df,aes(x=TIME,y=PD_lag,color=ID))+
  geom_line() +
  xlab("Time after dose (h)")+
  ylab("Effect")+
  theme_bw()+ 
  theme(legend.position = "none") 

Now, if we want to identify if there is a delay in the response after dosing, the most simple solution would be to replace the x=TIME with x=CP. This will plot the concentration-effect relationship, right? Unfortunately, this is wrong as can be observed below:

ggplot(df,aes(x=CP,y=PD_lag,color=ID))+
  geom_line() +
  xlab("Concentration")+
  ylab("Effect")+
  theme_bw()+ 
  theme(legend.position = "none")  

It does not give any error or warning when creating this plot. However, what happens is that the dots are being connected within an individual at increasing concentrations. The secret function we need to use here is geom_path(). With geom_path(), the data is plotted on a row-by-row basis and it connects the observations in chronological order (
https://ggplot2.tidyverse.org/reference/geom_path.html ).

# Effect with delay
p4<- ggplot(df,aes(x=CP,y=PD_lag,color=ID))+
  geom_path()+
  xlab("Concentration")+
  ylab("Effect")+
  theme_bw()+ 
  theme(legend.position = "none") 


# Direct effect
p5<-ggplot(df,aes(x=CP,y=PD_direct,color=ID))+
  geom_path()+
  xlab("Concentration")+
  ylab("Effect")+
  theme_bw()+ 
  theme(legend.position = "none") 


ggsave("plot_correct_hysteresis.png",plot=plot_grid(p4,p5),width=10)

 
Left = indirect effect, right = direct effect

This already provides us with much more information! However, it becomes quite messy with overlapping individual lines and summarizing the data per sampling time can help.

Summarize the data for a hysteresis figure

The difficulty in plotting the mean and standard deviation in a hysteresis figure is that there is variability in both directions at each timepoint. We have variability in the y axis, the effect size, and in the x-axis, the plasma concentrations.

We can calculate these values at each timepoint as follows:

########################### Calculate mean and sd in y and x direction
PD_sum <- df %>%
  group_by(TIME) %>%
  summarise(PD_mean = mean(PD_lag,na.rm=T),
            PD_sd=sd(PD_lag,na.rm=T),
            C_mean = mean(CP,na.rm=T),
            C_sd=sd(CP,na.rm=T)) %>% as.data.frame()
 

We can then plot the mean at each timepoint:

plot_pd_sum <- ggplot(PD_sum ,aes(x=C_mean,y=PD_mean))+
  geom_path(size=1.5)+
  theme_bw()+
  xlab("Concentration")+
  ylab("Effect") 

In which we can add the errorbars. Note: we need the geom_errorbarh() function to add the horizontal bars:

## Add error bars
plot_pd_sum_error <- plot_pd_sum+
  geom_errorbar(aes(ymin=PD_mean-PD_sd, ymax=PD_mean+PD_sd)) +
  geom_errorbarh(aes(xmin=C_mean-C_sd, xmax=C_mean+C_sd)) 

Now, if we do the same for the direct effect PD response we get the following figure:

You can observe the clear difference in hysteresis and it gives a clear overview compared to the individual lines to better inform decision making.

Conclusion

Creating hysteresis figures can be a piece of cake, once you have the right code. I hope this post provides you with the code needed to quickly create your own PK/PD hysteresis figures in R.

COMMENT

Any suggestions or typo’s? Leave a comment or contact me at info@pmxsolutions.com!

Leave a Reply

Your email address will not be published. Required fields are marked *