PK with deSolve

pkpd.Rmd

---
title: "Pharmacokinetics Pharmacodynamics with `DeSolve`"
author: "INSA Lyon - 3BIM EDO Modelling"
date: '2019-01-28'
output: slidy_presentation
---

```{r, message=FALSE, warning=FALSE, include=FALSE}
library(DiagrammeR)
library(deSolve) # need library `deSolve`
```

## PK/PD

Pharmacokinetics (PK)
:   The branch of pharmacology dedicated to determining the fate of substances administered to a living organism

Pharmacodynamics (PD)
:   The branch of pharmacology dedicated to study of the biochemical and physiologic effects of drugs, and in particular, pharmaceutical drugs.

The objective is to describe the *exposure-response* relationship.

## PK Models

How the body affects the a chemical (xenobiotic) after administration, through *absorption*, *distribution*, *metabolic changes*, and *excretion of metabolites*. *Multi-compartmental models* are the most common models. 

*ADME* or (*LADME*) scheme: Liberation, Absorption, Distribution, Metabolism, Excretion

Assumptions

- Immediate distribution within the compartment (drug concentration constant within a compartment)
- elimination and distribution rates are linear (only proportional to drug concentrations)



## Monocomparmental Models

*Monocompartmental* and *two compartmental* models are the most frequently used.

- Monocompartmental model: single physiological compartment, usually *plasma*, single drug (no metabolites). 

$$C(t) = C_0 \exp(- k_e t).$$

$C$ ($\mu$g/L) is the drug concentration in plasma, $C_0$ is the initial concentration. $C_0$ is related to the *dose* $D$ ($\mu$g) and the *volume of distribution* $V_d$ (L):

$$C_0 = \frac{D}{V_d}.$$
$k_e$ is the drug *elimination* rate constant.  

## Nomenclature {.smaller}

Term                    Description               Symbol        Equation          Example
------                  ------------              -------       ---------         --------
Dose                    Amount of drug            $D$           input par         250mg
                        administered 
Volume of               Apparent volume           $V_d$         $D/C_0$           5.0L
distribution            in which the drug 
                        is distributed
Concentration           Initial concentration     $C_0$         $D/V_d$           0.2mg/L
                        of the drug               
Elimination             Rate at which drug        $k_e$                           0.12h${}^{-1}$
rate constant           is removed from the
                        body
Area under              Integral of               $AUC_{0,t}$   $\int_0^tC(t)dt$  mg/L h
the curve (AUC)         drug concentration
                        over time 
Clearance               Volume of plasma          $CL$          $V_d k_e$         L/h
                        cleared of the drug
                        per unit time 

## Monocomparmental Models - Single Dose

The model can be expressed as an initial value problem: solve 

$$C' = - k_e C, \; C(0) = C_0.$$

```{r, message=FALSE, warning=FALSE}
monocompartment <- function(t, var, pars) {
  with(as.list(c(pars, var)), {
    dC <- - ke * C
    list(dC) }) }
D <- 500 # Dose (in mg)
Vd <- 5.0 # Volume of distribution (in L)
C0 <- c(C = D/Vd) # Initial concentration (in mg/L)
times <- seq(0, 100, 0.1) # (in hours)
pars <- c(ke = 0.04) # elimination rate (in h^-1)
sol <- ode(times = times, y = C0, 
           func = monocompartment, parms = pars)
```

## Monocomparmental Models - Single Dose

```{r, fig.width=6, fig.height=5}
plot(sol, ylab="concentration (mg/L)", ylim = c(0,D/Vd))
```

## Monocomparmental Models - Multiple Doses

- The elimination rate constant can include several elimination pathways: $k_e = k_{renal} + k_{met} + ...$.
- Multiple doses at dosing interval $\tau$ can be solved with the `events` function

```{r}
# 500mg per dose, three times a day 
# (08:00, 12:00, 18:00) for 2 days
dosingSchedule <- data.frame(var = rep(c("C"), 6),
  time = rep(c(8,12,18),2) + rep(c(0,24), each = 3),
  value = rep(c(D/Vd), 6),
  method = rep(c("add"), 6))
C0 <- c(C = 0) # Initial concentration = 0 mg/L
sol <- ode(times = times, y = C0, 
  func = monocompartment, parms = pars, 
  event = list(data = dosingSchedule))
```

## Monocomparmental Models - Multiple Doses

```{r, fig.width=6, fig.height=4}
plot(sol, ylab="concentration (mg/L)")
with(dosingSchedule,points(time, 0*time, pch = 19))
```

## Two compartmental Models - Formula

- Monocompartmental models do not take into account exchanges between plasma and peripheral organs 
- Two compartmental model: two physiological compartments, *central* (usually *plasma* or well perfused tissue) and the peripheral tissue (muscles, fat, brain, ...). The concentration at time $t$ in the central compartment is 

$$C(t) = a \exp(- \alpha t) + b \exp(- \beta t), \; \alpha > \beta.$$


- Alpha-phase: short-term kinetics during which the drug is distributed to the peripheral compartment
- Beta-phase: long-term kinetics during which the drug is eliminated

## Two compartmental Models - Alpha and Beta phases {.smaller}

```{r fig.width=5, fig.height=3}
pars <- c(a=80.0, b=20.0, alpha=1.2, beta=0.02)
t <- seq(0,40,by=0.1)
with(as.list(pars), 
  plot(t, a*exp(-alpha*t) + b*exp(-beta*t),
    type = 'l', log = 'y', lwd = 2, 
    xlab = 'time (h)', ylab = 'central compartment'))
```

## Two compartmental Models - Alpha and Beta phases {.smaller}

```{r, echo=FALSE, fig.width=6, fig.height=5}
with(as.list(pars), 
  plot(t, a*exp(-alpha*t) + b*exp(-beta*t),
    type = 'l', log = 'y', lwd = 2, 
    xlab = 'time (h)', ylab = 'central compartment'))
with(as.list(pars), lines(t, (a+b)*exp(-alpha*t), col = 'blue'))
with(as.list(pars), lines(t, b*exp(-beta*t), col = 'red'))
```

## Two compartmental Models - ODE formulation

```{r echo=FALSE}
DiagrammeR::grViz("
digraph rmarkdown  {
graph [splines=line]
Administration -> Central [label=\"D\"];
Central -> Elimination [label=\"k10\"];
Central -> Peripheral [arrowhead=rvee label=\"k12\" color=\"black:white:black\" labeldistance=20];
Peripheral -> Central [arrowhead=rvee label=\"k21\" color=\"black:white:black\" labeldistance=20];
}                  
", height = 200)
```

The pair of ODEs is

\begin{align*}
C_c' & = -k_{12} C_c + k_{21} C_p - k_{10} C_c \\
C_p' & =  k_{12} C_c - k_{21} C_p \\
\end{align*}

Parameters $a,b,\alpha,\beta$ and $D,k_{12},k_{21},k_{10}$ are related.

## Two compartmental Models - ODE formulation

In the two compartment model, the *volume of distribution* changes as the drug gets distributed into the 
peripheral compartment. 
We denote $V_c,V_p$ the volumes of distribution
of the central and peripheral compartments. 

```{r, message=FALSE, warning=FALSE}
twocompartment <- function(t, var, pars) {
  with(as.list(c(pars, var)), {
    dCc <- -k12*Central + k21*r*Peripheral - k10*Central
    dCp <-  k12/r*Central - k21*Peripheral
    list(c(dCc,dCp)) }) }
D <- 500 # Dose (in mg)
Vc <- 5.0 # Volume of distribution in Central (in L)
Vp <- 10.0 # Volume of distribution in Peripheral (in L)
C0 <- c(Central = D/Vc, Peripheral = 0) # Initial concentration (in mg/L)
times <- seq(0, 40, 0.1) # (in hours)
pars <- c(k12 = 1.2, k21 = 1.5, k10 = 0.04, r = Vp/Vc) 
sol <- ode(times = times, y = C0, 
           func = twocompartment, parms = pars)
```

## Two compartmental Models - ODE formulation

```{r, fig.width=6, fig.height=5}
plot(sol, which = 'Central',
    type = 'l', log = 'y', lwd = 2, 
    xlab = 'time (h)', ylab = 'central compartment')
```

This is cool. I am trying to set up an "interface model for dosage adjustment connects hematotoxicity to pharmacokinetics". In this model there is a feedback loop and the influx must not only be amplified but also be delayed. Is there any good starting point to implement such models Image and equations Shown below.

twocompartment <- function(t, var, pars) {
  with(as.list(c(pars, var)), {
    dCc <- -k12*Central + k21*r*Peripheral - k10*Central
    dCp <-  k12/r*Central - k21*Peripheral
     dCy <- -alpha * exp(-beta * cumexp) * cumexp + (Central - gamma) *H(Central - gamma) 
    dCs <- ks*kw*((neut/mature)^-0.2)*(1-exp(-u*cumexp))- ks * neut
    dCw <- -kw * mature + (t-k) * neut
    
    list(c(dCc,dCp,dCy,dCs,dCw)) }) }
D <- 500 # Dose (in mg)
Vc <- 5.0 # Volume of distribution in Central (in L)
Vp <- 10.0 # Volume of distribution in Peripheral (in L)
C0 <- c(Central = D/Vc, Peripheral = 0,cumexp=0,neut=0,mature=3) # Initial concentration (in mg/L)


times <- seq(0, 40, 0.1) # (in hours)
pars <- c(k12 = 1.2, k21 = 1.5, k10 = 0.04, r = Vp/Vc,Rs = 1, d = 0.1, k = 1.65, kw = 1, ks = 0.5, theta = 1, W = 100, K = 100,
          alpha = 0.85, beta = 1.5*10-4, gamma = 0.2,u=5.6*10-4,s=s0) 
sol <- ode(times = times, y = C0, 
           func = twocompartment, parms = pars)

my model does not simulate for some reason.
below is the source for this. The only difference is 2 vs 3 compartment

https://link.springer.com/article/10.1007/s10928-008-9106-4#Sec2

The parameter s=s0 is unused, and s0 is undefined, consider removing s from the list of parameters pars.
The function H is not defined. If this is the heaviside function, you can use the logical comparison ( (Central - gamma) >= 0 ), which evaluates to 1 if Central is larger than gamma.
You will also have a problem with the term with a negative exponent (neut/mature)^-0.2; when neut is 0.0, this will evaluate to Inifinity.

There are two ways you can delay the feedback loop. If your feedback depends on a variable z, you can

1. Introdude a linear chain:

    dx1/dt = beta*(z - x1)
    dx2/dt = beta*(x1 - x2)
    dx3/dt = beta*(x2 - x3)

The variable x3 will be a lagged version of z. This method increase the number of variables in your system.

2. Use an explicit delay, by evaluating z at t-tau: z(t-tau). To solve differential equations with delays, you will need to use the delay differential equation solver dede in the library deSolve.