Automatic model diagrams

A compartment diagram is often the fastest way to explain a model in a report or to check that the equations say what you meant. nlmixr2plot can draw one directly from a model’s differential equations:

A first diagram

Any model function, rxode2 model, rxode2 user interface object or fitted nlmixr2 object can be diagrammed. Here is a two-compartment model with first-order absorption:

two.cmt <- function() {
  ini({
    tka <- log(1.5)
    tcl <- log(3)
    tv <- log(20)
    tq <- log(2)
    tvp <- log(40)
    add.sd <- 0.2
  })
  model({
    ka <- exp(tka)
    cl <- exp(tcl)
    v <- exp(tv)
    q <- exp(tq)
    vp <- exp(tvp)
    d/dt(depot) <- -ka * depot
    d/dt(central) <- ka * depot - cl / v * central - q / v * central +
      q / vp * periph
    d/dt(periph) <- q / v * central - q / vp * periph
    cp <- central / v
    cp ~ add(add.sd)
  })
}

modelDiagram(two.cmt, engine = "ggplot2")
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments

plot() of any rxode2 user interface object (or compiled rxode2 model) draws the same diagram, so this is equivalent:

plot(rxode2::rxode2(two.cmt), engine = "ggplot2")

The dosing compartment (depot, drawn with a heavy border) is on top and feeds the central compartment. The peripheral compartment, which exchanges mass with central in both directions, is to the left, and elimination leaves the central compartment downwards.

Drawing engines

The same graph can be drawn several ways, chosen with engine:

"DiagrammeR" draws an interactive Graphviz widget with the ‘DiagrammeR’ package. This is the default when ‘DiagrammeR’ is installed.
"ggplot2" returns a ggplot object that can be themed, saved with ggplot2::ggsave() or combined with other plots. It needs no extra packages.
"dot" returns the Graphviz DOT source as a string, to edit or to render with any Graphviz tool.

modelDiagram(two.cmt, engine = "DiagrammeR")
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments

cat(modelDiagram(two.cmt, engine = "dot"))
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> digraph model {
#>   graph [layout = neato, splines = true, outputorder = edgesfirst, forcelabels = true];
#>   node [shape = box, style = "rounded,filled", fontname = Helvetica];
#>   edge [fontname = Helvetica, fontsize = 10];
#>   "depot" [pos = "0,1.1!", fillcolor = "#F2C57C", penwidth = 2];
#>   "central" [pos = "0,0!", fillcolor = "#7FB3D5"];
#>   "periph" [pos = "-1.67,0!", fillcolor = "#A9CCE3"];
#>   "depot" -> "central";
#>   "central" -> "periph" [dir = both];
#>   ".elimination4" [shape = point, style = invis, width = 0.01, pos = "0,-0.77!"];
#>   "central" -> ".elimination4";
#> }

The default engine can be set once per session:

options(nlmixr2plot.diagram.engine = "ggplot2")

Arrows can be labeled with the model terms that drive them:

modelDiagram(two.cmt, engine = "ggplot2", labels = TRUE)
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments

The graph behind the diagram

modelGraph() returns the parsed graph, which is useful for checking how the equations were interpreted (or for drawing it some other way):

g <- modelGraph(two.cmt)
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
g
#> nlmixr2 model graph
#> 
#> compartments:
#>     name       role dosing
#>    depot     dosing   TRUE
#>  central    central  FALSE
#>   periph peripheral  FALSE
#> 
#> flows:
#>     from       to        type sign          label
#>    depot  central    transfer    1     ka * depot
#>  central   periph    transfer    1  q/v * central
#>   periph  central    transfer    1  q/vp * periph
#>  central (output) elimination   -1 cl/v * central

The nodes data frame gives each compartment’s role and layout position, and the edges data frame lists every flow. A graph can be passed straight to modelDiagram() or plot():

plot(g, engine = "ggplot2")

How the equations are read

Each d/dt() equation is split into signed additive terms (products are distributed over sums), and then:

Mass transfer: a term subtracted from one compartment and added, identically, to another moves mass between them. The order of factors does not matter (ka*depot matches depot*ka), and first-order, zero-order and enzyme-driven rates are all recognized.
Elimination: a remaining loss that contains the compartment’s own amount (like cl/v*central or Michaelis-Menten vmax*C/(km + C)), or that depends on no compartment at all (a zero-order loss).
Input: a remaining gain that depends on no other compartment: a zero-order kin, or self-dependent growth like kg*A in d/dt(A).
Interaction: a term that depends on another compartment without moving mass, like an effect compartment or a pharmacodynamic stimulation or inhibition. It is drawn dashed. Its direction comes from the equation: inhibition when the term decreases as the driving compartment increases (including 1 - emax*C/(ec50 + C), kin/(1 + C) and exp(-k*C) forms), stimulation when it increases, and “modulation” when the direction cannot be determined.

Intermediate variables (like cp <- central/v) are followed, if/else blocks and ifelse() keep their conditions, and residual error lines are ignored.

Pharmacodynamic models

Compartments that interact with the pharmacokinetic model without mass transfer go to the right, with their own inputs above and outputs below. Here is an indirect response (turnover) model where the drug inhibits the production of the response:

pk.turnover <- function() {
  ini({
    tktr <- log(1)
    tka <- log(1)
    tcl <- log(0.1)
    tv <- log(10)
    poplogit <- 2
    tec50 <- log(0.5)
    tkout <- log(0.05)
    te0 <- log(100)
    prop.err <- 0.1
    pdadd.err <- 10
  })
  model({
    ktr <- exp(tktr)
    ka <- exp(tka)
    cl <- exp(tcl)
    v <- exp(tv)
    emax <- expit(poplogit)
    ec50 <- exp(tec50)
    kout <- exp(tkout)
    e0 <- exp(te0)
    DCP <- center / v
    PD <- 1 - emax * DCP / (ec50 + DCP)
    effect(0) <- e0
    kin <- e0 * kout
    d/dt(depot) <- -ktr * depot
    d/dt(gut) <- ktr * depot - ka * gut
    d/dt(center) <- ka * gut - cl / v * center
    d/dt(effect) <- kin * PD - kout * effect
    cp <- center / v
    cp ~ prop(prop.err)
    effect ~ add(pdadd.err)
  })
}
modelDiagram(pk.turnover, engine = "ggplot2")
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments

The transit compartment gut stacks between the dosing compartment and the central compartment. The response has a zero-order input (kin) above it, an elimination below it, and a dotted arrow from center: the drug inhibits the response, because PD decreases as the concentration increases. With ‘DiagrammeR’ the inhibition is drawn with a “tee” arrow head:

modelDiagram(pk.turnover, engine = "DiagrammeR")
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments

A more involved example, with a metabolite, two peripheral compartments, an effect compartment and a response driven by the effect compartment:

pkpd <- rxode2::rxode2({
  C2 <- centr / V2
  C3 <- peri / V3
  C4 <- peri2 / V4
  d/dt(depot) <- -KA * depot
  d/dt(centr) <- KA * depot - CL * C2 - Q * C2 + Q * C3 - Q2 * C2 + Q2 * C4 -
    kmet * centr
  d/dt(peri) <- Q * C2 - Q * C3
  d/dt(peri2) <- Q2 * C2 - Q2 * C4
  d/dt(met) <- kmet * centr - kelm * met
  d/dt(ce) <- ke0 * (C2 - ce)
  d/dt(resp) <- kin - kout * (1 - ce / (ec50 + ce)) * resp
})
modelDiagram(pkpd, engine = "ggplot2")

Target-mediated drug disposition

Binding and unbinding are recognized as mass transfer from both binding partners into the complex (and back):

tmdd <- rxode2::rxode2({
  d/dt(central) <- -kel * central - kon * central * target + koff * complex
  d/dt(target) <- ksyn - kdeg * target - kon * central * target +
    koff * complex
  d/dt(complex) <- kon * central * target - koff * complex - kint * complex
})
modelGraph(tmdd, dosing = "central")
#> nlmixr2 model graph
#> 
#> compartments:
#>     name       role dosing
#>  central    central   TRUE
#>   target peripheral  FALSE
#>  complex peripheral  FALSE
#> 
#> flows:
#>     from       to        type sign                  label
#>  central  complex    transfer    1 kon * central * target
#>   target  complex    transfer    1 kon * central * target
#>  complex  central    transfer    1         koff * complex
#>  complex   target    transfer    1         koff * complex
#>  central (output) elimination   -1          kel * central
#>  (input)   target       input    1                   ksyn
#>   target (output) elimination   -1          kdeg * target
#>  complex (output) elimination   -1         kint * complex
modelDiagram(tmdd, dosing = "central", engine = "ggplot2")

Dosing compartments

Dosing compartments are detected from the dosing records of the data: for a fitted model this is the data it was fit to, and for other models it can be supplied with data. Both numeric and named cmt values are understood:

d <- data.frame(
  id = 1, time = c(0, 0, 1, 2),
  amt = c(100, 50, 0, 0), evid = c(1, 1, 0, 0),
  cmt = c("depot", "central", "central", "central"), dv = 0
)
modelGraph(two.cmt, data = d)$nodes
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#>      name       role dosing  x y annotation
#> 1   depot     dosing   TRUE  0 1           
#> 2 central    central   TRUE  0 0           
#> 3  periph peripheral  FALSE -1 0

With no data, the first compartment (rxode2’s default dosing compartment) is assumed. The dosing compartments can also be given directly:

modelDiagram(two.cmt, dosing = "central", engine = "ggplot2")
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments

For a fit, modelDiagram(fit) uses the fitted model and data:

fit <- nlmixr2(two.cmt, nlmixr2data::theo_sd, est = "focei")
modelDiagram(fit)

Dosing properties and delays

Dosing properties set in the model (lag()/alag(), f()/F(), rate() and dur()) do not change the flows, so they are shown as an annotation next to their compartment (compartments without them are left blank). Initial conditions like central(0) <- 0 do not change the diagram.

pk.lag <- rxode2::rxode2({
  d/dt(depot) <- -ka * depot
  alag(depot) <- tlag
  f(depot) <- fbio
  d/dt(central) <- ka * depot - cl / v * central
  dur(central) <- d1
})
modelGraph(pk.lag)
#> nlmixr2 model graph
#> 
#> compartments:
#>     name    role dosing           annotation
#>    depot  dosing   TRUE lag = tlag; F = fbio
#>  central central  FALSE             dur = d1
#> 
#> flows:
#>     from       to        type sign          label
#>    depot  central    transfer    1     ka * depot
#>  central (output) elimination   -1 cl/v * central
modelDiagram(pk.lag, engine = "ggplot2")

A delay() keeps the meaning of what it delays: ka*delay(depot, tlag) in the destination still matches -ka*depot in the source (a delayed transfer), and a delayed concentration driving an effect keeps the direction of the concentration.

`linCmt()` models

Models written with linCmt() are converted to their ODE form with rxode2::linToOde() before they are diagrammed (this requires a version of ‘rxode2’ that provides it):

one.cmt <- function() {
  ini({
    tka <- 0.45
    tcl <- 1
    tv <- 3.45
    add.sd <- 0.7
  })
  model({
    ka <- exp(tka)
    cl <- exp(tcl)
    v <- exp(tv)
    linCmt() ~ add(add.sd)
  })
}
modelDiagram(one.cmt, engine = "ggplot2")
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments

Limitations

The diagram is only as good as the parsing of the equations, so a few conventions matter:

Mass transfer is only detected when the same term is subtracted from the source and added to the destination. A transfer scaled in only one of the equations (for example a volume or stoichiometric conversion) is drawn as an elimination plus an interaction.
A gain driven only by another compartment (like ke0*cp into an effect compartment) is represented by the dashed interaction arrow alone.
The layout is automatic; for publication-quality tweaks, start from the DOT source (engine = "dot") or the ggplot object and adjust it.