# Personalizing drug doses in the cloud using R

## > Sys.getenv(country = ‘USA’)

If medications were prescribed, monitored and taken properly,
we wouldn’t face this cost, and patients would be healthier. Watanabe et al, 2018 (doi.org/10.1177/10600280187651)

## > demo(‘rx.studio’)

youtu.be/MZUrMFbd9TA

## > vignette(topic = ‘PK/PD models’)

Source: Mortensen et al (2008): Introduction to PK/PD modelling.

## > vignette(topic = ‘PK/PD models’)

Source: Mortensen et al (2008): Introduction to PK/PD modelling.

## > vignette(topic = ‘PK/PD models’)

Source: Mortensen et al (2008): Introduction to PK/PD modelling.

## > vignette(topic = ‘PK/PD models’)

Source: Mortensen et al (2008): Introduction to PK/PD modelling.

## > vignette(topic = ‘PK/PD models’)

Source: Spherical cow

## > vignette(topic = ‘PK/PD models’)

$C=\frac{A}{V}$

• $$C$$ drug concentration
• $$A$$ drug amount
• $$V$$ volume of distribution

Example: 500 mg Panadol ($$Vd = 0.9L/kg$$) administered for a 70 kg patient

$C=\frac{500mg}{70kg * 0.9L/kg}=\frac{500mg}{63L}=7.9 mg/L$

## > vignette(topic = ‘PK/PD models’)

$C_{oral}(t)=\frac{A_{oral}(t)}{V}=\frac{K_aFA_0}{V(K_a-K)}(exp(-K \cdot t) - exp(-K_a \cdot t))$

• $$C$$ drug concentration
• $$A$$ drug amount
• $$V$$ volume of distribution
• $$K_a$$ absorption constant
• $$K$$ elimination rate
• $$F$$ bioavailability

## > vignette(topic = ‘PK/PD models’)

$C_{oral}(t)=\frac{A_{oral}(t)}{V}=\frac{K_aFA_0}{V(K_a-K)}(exp(-K \cdot t) - exp(-K_a \cdot t))$

#' Concentration at a time computed using a one-compartment model (oral dose)
#' @param t time (hours)
#' @param dose dose amount (mg)
#' @param v volume of distribution (l)
#' @param k elimination rate constant (h^-1)
#' @param ka absorption rate constant (h^-1)
#' @param f bioavailability
#' @return numeric
#' @export
ct <- function(t, dose, v, k, ka, f) {
(ka * f * dose) / (v * (ka - k)) * (exp(-k * t) - exp(-ka * t))
}

## > example(topic = ‘paracetamol’)

Source: Rawlins et al. (1977): Paracetamol (simplified)

## > example(topic = ‘paracetamol’)

Source: Rawlins et al. (1977): Paracetamol (simplified)

## > example(topic = ‘paracetamol’)

#' Concentration at a time computed using a one-compartment model (oral dose)
#' @param t time (hours)
#' @param dose dose amount (mg)
#' @param v volume of distribution (l)
#' @param k elimination rate constant (h^-1)
#' @param ka absorption rate constant (h^-1)
#' @param f bioavailability
#' @return numeric
#' @export
ct <- function(t, dose, v, k, ka, f) {
(ka * f * dose) / (v * (ka - k)) * (exp(-k * t) - exp(-ka * t))
}
ctp <- purrr::partial(ct, v = 42, k = 0.28, ka = 1.8, f = 0.89) # 70 kgs adult
> ctp(t = 1, dose = 1000)
[1] 14.81762

## > example(topic = ‘paracetamol’)

#' Concentration at a time computed using a one-compartment model (oral dose)
#' @param t time (hours)
#' @param dose dose amount (mg)
#' @param v volume of distribution (l)
#' @param k elimination rate constant (h^-1)
#' @param ka absorption rate constant (h^-1)
#' @param f bioavailability
#' @return numeric
#' @export
ct <- function(t, dose, v, k, ka, f) {
(ka * f * dose) / (v * (ka - k)) * (exp(-k * t) - exp(-ka * t))
}
ctp <- purrr::partial(ct, v = 42, k = 0.28, ka = 1.8, f = 0.89) # 70 kgs adult
> ctp(t = 1, dose = 1000)
[1] 14.81762
> ctp(t = 2, dose = 1000)
[1] 13.64825

## > example(topic = ‘paracetamol’)

#' Concentration at a time computed using a one-compartment model (oral dose)
#' @param t time (hours)
#' @param dose dose amount (mg)
#' @param v volume of distribution (l)
#' @param k elimination rate constant (h^-1)
#' @param ka absorption rate constant (h^-1)
#' @param f bioavailability
#' @return numeric
#' @export
ct <- function(t, dose, v, k, ka, f) {
(ka * f * dose) / (v * (ka - k)) * (exp(-k * t) - exp(-ka * t))
}
ctp <- purrr::partial(ct, v = 42, k = 0.28, ka = 1.8, f = 0.89) # 70 kgs adult
> ctp(t = 1, dose = 1000)
[1] 14.81762
> ctp(t = 2, dose = 1000)
[1] 13.64825
> ctp(t = 6, dose = 1000)
[1] 4.676354

## > example(topic = ‘paracetamol’)

library(data.table); library(ggplot2)
conc <- data.table(h = seq(0, 24, by = 0.1))
conc[, c := ctp(h, 1000)]
ggplot(conc, aes(h, c)) + geom_line()

## > do.call(ctp, weights)

ctp <- function(t, dose, weight, k = 0.28, ka = 1.8, f = 0.89) {
ct(t, dose, v = weight * 0.6, k, ka, f)
}
conc[, c1 := ctp(h, 1000, weight =  40)]
conc[, c2 := ctp(h, 1000, weight =  60)]
conc[, c3 := ctp(h, 1000, weight =  80)]
conc[, c4 := ctp(h, 1000, weight = 100)]
ggplot(melt(conc, id.vars = 'h'), aes(h, value, color = variable)) +
geom_line() +
scale_color_discrete(name = 'Weight',
labels = c('40 kg', '60 kg', '80 kg', '100kg')) +
ylab('Blood concentration forecast') +
theme(legend.position = 'top') +
geom_hline(yintercept = 10, color = 'black', linetype = 2, size = 1.25) +
geom_hline(yintercept = 20, color = 'black', linetype = 2, size = 1.25)

## > do.call(ctp, family)

An anonymized family:

• GD (96 kg male): 1 pill of Panadol (500 mg)
• HD (64 kg female): 1 pill of Panadol (500 mg)
• BD (41 kg male): 1 pill of Panadol (500 mg)
• BD (23 kg female): half pill of Panadol (250 mg)
• BD (13 kg female): 7-85ml Panadol baby (120 mg)

Therapeutic goals for paracetamol:

• sub-therapeutic range: < 10 mg/L concentrations
• therapeutic range: ≥ 10 and ≤ 20 mg/L concentrations
• toxic: > 75 mg/L concentrations

## > do.call(ctp, family)

conc[, c1 := ctp(h, 500, weight = 96)]
conc[, c2 := ctp(h, 500, weight = 64)]
conc[, c3 := ctp(h, 500, weight = 41)]
conc[, c4 := ctp(h, 250, weight = 23)]
conc[, c5 := ctp(h, 120, weight = 13)]
ggplot(melt(conc, id.vars = 'h'), aes(h, value, color = variable)) + geom_line() +
scale_color_discrete(name = 'Weight', labels = c('DG', 'DH', 'DB', 'DB', 'DB')) +
theme(legend.position = 'top') +
geom_hline(yintercept = 10, color = 'black', linetype = 2, size = 1.25)

## > do.call(ctp, family, ndoses = 4)

$C_{MD}(t)=\sum_{n=0}^{N-1} C_{oral}(t-nτ)$

• $$N$$ number of doses
• $$τ$$ dosing interval

Source: Mortensen et al (2008): Introduction to PK/PD modelling.

ctpm <- function(t, doses, interval, dose, weight) {
sum(sapply(0:(doses-1), function(n) ctp(t - n*interval, dose, weight)))
}

## > do.call(ctp, family, ndoses = 4)

ctpm <- function(t, doses, interval, dose, weight) {
sum(sapply(0:(doses-1), function(n) ctp(t - n*interval, dose, weight)))
}

conc <- data.table(h = seq(0, 24, by = 0.1))
conc[, doses := pmin(h %/% 4, 3) + 1]
conc[, c1 := ctpm(h, doses, 4, 500, weight =  96), by = .(h, doses)]
conc[, c2 := ctpm(h, doses, 4, 500, weight =  64), by = .(h, doses)]
conc[, c3 := ctpm(h, doses, 4, 500, weight =  41), by = .(h, doses)]
conc[, c4 := ctpm(h, doses, 4, 250, weight =  23), by = .(h, doses)]
conc[, c5 := ctpm(h, doses, 4, 120, weight =  13), by = .(h, doses)]

ggplot(melt(conc[, -'doses'], id.vars = 'h'), aes(h, value, color = variable)) +
geom_line() +
scale_color_discrete(name = 'Weight', labels = c('DG', 'DH', 'DB', 'DB', 'DB')) +
theme(legend.position = 'top') +
geom_hline(yintercept = 10, color = 'black', linetype = 2, size = 1.25) +
geom_hline(yintercept = 20, color = 'black', linetype = 2, size = 1.25)

## > sd(do.call(ctp, rlnorm(2000)))

Source: Rawlins et al. (1977): Paracetamol (simplified)

## > sd(do.call(ctp, rlnorm(2000)))

weight <- 70
meanlog <- log((weight * 0.6)^2 / sqrt(0.07^2 + (weight * 0.6)^2))
sdlog <- sqrt(log(1 + (0.07^2 / (weight * 0.6)^2)))
ggplot(data.frame(x = rlnorm(n = 2000L, meanlog, sdlog))) +
geom_histogram(aes(x)) +
ggtitle('Volume of distribution for 70 kg') + xlab('')

## > sd(do.call(ctp, rlnorm(2000)))

for (i in 1:250) {
meanlog <- log((weight * 0.6)^2 / sqrt(0.07^2 + (weight * 0.6)^2))
sdlog <- sqrt(log(1 + (0.07^2 / (weight * 0.6)^2)))
conc <- copy(conc)[, c := ct(h, 1000, v = rlnorm(n = 1L, meanlog, sdlog),
k = 0.28, ka = 1.8, f = 0.89)]
(G <- G + geom_line(data = conc, color = 'gray', alpha = 0.1))
}
G + geom_line(color = 'black') + ggtitle('250 simulations')

## > sd(do.call(ctp, rlnorm(2000)))

simdata <- rbindlist(lapply(1:2000, function(i) {
meanlog <- log((weight * 0.6)^2 / sqrt(0.07^2 + (weight * 0.6)^2))
sdlog <- sqrt(log(1 + (0.07^2 / (weight * 0.6)^2))) * 5
as.data.frame(matrix(
ct(seq(0, 24, by = 0.1), 1000, v = rlnorm(n = 1L, meanlog, sdlog),
k = 0.28, ka = 1.8, f = 0.89),
nrow = 1))
}))

simagg <- data.frame(
h = seq(0, 24, by = 0.1),
min = apply(simdata, 2, FUN = min),
mean = apply(simdata, 2, FUN = mean),
max = apply(simdata, 2, FUN = max))

ggplot(simagg, aes(h)) +
geom_ribbon(aes(ymin=min, ymax=max), fill = 'gray') +
geom_line(aes(y = mean)) +
ylab('Blood concentration forecast') +
theme_bw() + theme(legend.position = 'top')

## > fit(ctp, data.frame(conc = 5))

library(FME)

concentrations <- data.table(x = 5, y = 5)
model_function <- function(params, x) {
ct(x, dose = 1000, v = params['v'], k = 0.28, ka = 1.8, f = 0.89)
}
cost_function <- function(params) {
out <- model_function(params, concentrations$x) concentrations$y - out
}

bay <- modFit(f = cost_function, p = c(v = weight * 0.6), method = 'Newton')
bayc <- copy(conc)[, c := ct(seq(0, 24, by = 0.1), 1000, v = bay$par[['v']], k = 0.28, ka = 1.8, f = 0.89)] ggplot(simagg, aes(h)) + geom_ribbon(aes(ymin=min, ymax=max), fill = 'gray') + geom_line(aes(y = mean)) + geom_point(aes(x=5, y=5), color = 'orange', size = 3) + geom_line(data = bayc, aes(y = c), color = 'orange') + ylab('Blood concentration forecast') + theme_bw() + theme(legend.position = 'top') ## > fit(ctp, data.frame(conc = 5)) > 0.60 * 70 [1] 42 > bay$par[['v']]
[1] 51.954

## > library(chatgpt)

placeholder as cannot finish with an image