Mass balance in R's deSolve Using Events

Question

I reproduced a one-compartment model describing drug injection from "Solving Differential Equations in R" (see code below).

What is unclear to me is how to check that mass balance is 0 in case of events. Any advice?

library(deSolve)
b <- 0.6 
yini <- c(blood = 0)

pharmaco2 <- function(t, blood, p) { 
  dblood <-- b * blood 
  list(dblood) }

injectevents <- data.frame(var = "blood", time = 0:20, value = 40, method = "add")

times <- seq(from = 0, to = 10, by = 1/24) 

out2 <- ode(func = pharmaco2, 
            times = times, 
            y = yini, 
            parms = NULL, 
            method = "impAdams", 
            events = list(data = injectevents))

plot(out2, lwd = 2, xlab="days")

Answer 1

Checking mass balance requires to close the system. This can be done by adding explicit state variables for source and sink, where source represents the stock of the drug and sink the accumulated amount that leaves the blood, either by decomposition or excretion. The source variable can be set to zero to check mass balance, or any positive value to represent the amount of drug available in stock.

Then add corresponding processes to the model, i.e. decomposition/excretion as state equation and removal from the stock during injection in the events. The model reads then as follows:

library(deSolve)
b <- 0.6
yini <- c(source=0, blood = 0, sink=0)

pharmaco2 <- function(t, yini, p) {
  blood   <-  yini[2]
  dsource <-  0
  dblood  <- -b * blood
  dsink   <- +b * blood
  list(c(dsource, dblood, dsink))
}

injectevents <- data.frame(var = rep(c("blood", "source"), 21),
                           time = rep(0:20, each = 2),
                           value =  rep(c(40, -40), 21),
                           method = "add")

times <- seq(from = 0, to = 10, by = 1/24)

out2 <- ode(func = pharmaco2,
            times = times,
            y = yini,
            parms = NULL,
            method = "lsoda",
            events = list(data = injectevents))

plot(out2, lwd = 2, xlab="days", mfrow=c(1,3))

summary(rowSums(out2[,-1])) # sum of all columns, except time

And the mass balance check can be done by adding the values of all columns, except time. If initial source was set to zero, then the balance should also be zero, up to numerical rounding errors:

      Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
-1.990e-13 -4.974e-14 -2.309e-14 -2.541e-14  2.132e-14  9.948e-14 

Note also that I use default lsoda as a solver. Solver adamsor bdf work as well, and also impAdams but I would not recommend the latter in this case.