Strategies for Visualizing Hysteresis in PK/PD Data

What is hysteresis ?

Hysteresis in pharmacokinetic/pharmacodynamic relationships is a phenomenon arising when there is a time dependency on the relationship between drug concentrations and measured response. Hysteresis can be caused by a multitude of different reasons, most typically manifesting as either a delay between measured drug concentrations and down stream pharmacodynamic responses or a wearing off (e.g. tolerance) of drug effect over time. This results in different pharmacodynamic responses at equivalent drug concentrations depending on when in time the variables were measured.

Identifying Hysteresis in Your Data

To identify hysteresis in your PK/PD dataset, start by generating a plot of time-matched concentration vs. response variables. The points in the plot should be connected in the sequence they were measured (i.e. sorted by time). To create this plot in ggplot, we can use the function geom_path instead of the more commonly used geom_line. The main difference between the two is that geom_path will connect the points in the order in which they appear in the dataset, while geom_line will connect them from left to right. This also means that your data needs to be sorted by the time variable for the plot to display correctly.

Generate the PKPD data

We're going to use mrgsolve to generate some PKPD data from an indirect response model (specifically inhibition of production).

Simulating the data

library(mrgsolve)
library(dplyr)

mod <- mread("irm1", modlib())
dosing <- ev(ID = 1, amt = 25, time = 0, dose = amt) +
          ev(ID = 2, amt = 50, time = 0, dose = amt) +
          ev(ID = 3, amt = 75, time = 0, dose = amt) +
          ev(ID = 4, amt = 100, time = 0, dose = amt)

res <- mod %>%
  data_set(dosing) %>%
  param(CL = 0.5, KA = 0.1, KOUT = 0.1) %>%
  obsonly() %>%
  carry.out(dose) %>%
  mrgsim(end = 168, tgrid = c(0, 1, 2, 3, 4, 6, 8, 10,
                              12, 14, 16, 18, 21, 24,
                              28, 32, 36, 44, 52, 60,
                              72, 84, 96, 120, 144, 168)) %>%
  select(time, CP, RESP, ID, dose)

Take a peek at the dataset structure

print(head(res,10))

Plot the data

Let's plot concentration vs response to visualize the hysteresis.

library(ggplot2)

theme_set(theme_light())

plot1 <- ggplot(res, aes(x = CP, y = RESP, group = ID)) +
 geom_path(aes(colour = as.factor(dose))) +
 geom_point(aes(colour = as.factor(dose))) +
 scale_color_viridis_d() +
 facet_wrap('dose', nrow = 1) +
 theme(legend.position = "top") +
 labs(x = "Drug Concentration", y = "Response (% of Baseline)", color = "Dose (mg)")


print(plot1)

We can see a clear hysteresis, although it is difficult to tell in which direction the loop goes. When looking at your observed data, it is important to clearly visualize the direction of the hysteresis loop, as it can give you insights into what types of PK/PD models may be appropriate to use. A loop that shows more pharmacodynamic effect at later timepoints (at similar concentrations) suggests a delay in drug effect. A loop that shows less effect at later timepoints suggests a wearing off of drug effect over time.

We'll be using two strategies to illustrate which direction the hysteresis loop is going. The first will be to add a time label to each point and the second will be to add arrows to the plot in the direction of the loop.

Adding time labels

plot2 <- ggplot(res, aes(x = CP, y = RESP, group = ID)) +
  geom_path(aes(colour = as.factor(dose))) +
  geom_point(aes(colour = as.factor(dose))) +
  geom_text(aes(label = paste(time, "hr"),x = CP +0.2, y = RESP +1), size = 3) +
  scale_color_viridis_d() +
  facet_wrap('dose', nrow = 1) +
  theme(legend.position = "top") +
  labs(x = "Drug Concentration", y = "Response (% of Baseline)", color = "Dose (mg)")

print(plot2)
 

Adding arrows

To add arrows, we'll use geom_segment with an arrow segment end. To use geom_segment, we'll need to define a segment end in the dataset, where the end of the segment for each point will be the x and y variables of the next timepoint. The lead() function from dplyr can help us with this.

res <- res  %>%
  mutate(xend = lead(CP), yend = lead(RESP))

Next, we'll use geom_segment to create the plot we want.

plot3 <- ggplot(res, aes(x = CP, y = RESP, group = ID)) +
 geom_segment(aes(xend = xend, yend = yend, colour = as.factor(dose)),
              arrow = arrow(length = unit(0.08, "npc"), type="closed")) +
 scale_color_viridis_d() +
 facet_wrap('dose', nrow = 1) +
 theme(legend.position = "top") +
 labs(x = "Drug Concentration", y = "Response (% of Baseline)", color = "Dose (mg)")

print(plot3)
 
 

Using animation to visualize the hysteresis

In addition to the static plots created earlier, animation can also be a useful strategy for visualizing hysteresis data. We'll continue to plot a 2-dimension plot of concentration vs. response, except now we'll incorporate time as the third dimension represented by frames of the animation. To create this figure we'll use the gganimate package.

library(gganimate)

animation <- ggplot(res, aes(x = CP, y = RESP, group = ID)) +
  geom_point(aes(color = as.factor(dose)), size = 3) +
  geom_line(aes(color = as.factor(dose)), size = 1) +
  geom_text(aes(label = paste(time, "hr"), y = RESP +1)) +
  scale_color_viridis_d() +
  transition_reveal(time) +
  theme(legend.position = "top") +
  labs(x = "Drug Concentration", y = "Response (% of Baseline)",
       color = "Dose (mg)")
print(animation)

Have another way you like to visualize hysteresis? Let us know in the comments below.