Model Representation

Once your model is defined and your analyses run, you might want to show it off. Transcribing a computational model to equations can be a tedious affair. Furthermore, it is easy to get something wrong. We, therefore, sought to automate that process.

Models created by the Pumas.jl model macros contain all the information we need to automatically generate LaTeX equations from them using Latexify.jl.

Let us first define a model so that we can demonstrate this:

using Pumas

model = @model begin
  @param  begin
    tvcl ∈ RealDomain(lower=0)
    tvvc ∈ RealDomain(lower=0)
    tvq ∈ RealDomain(lower=0)
    tvvp ∈ RealDomain(lower=0)
    Ω_pk ∈ PDiagDomain(4)
    σ_prop_pk ∈ RealDomain(lower=0)
    # PD parameters
    tvturn ∈ RealDomain(lower=0)
    tvebase ∈ RealDomain(lower=0)
    tvec_50 ∈ RealDomain(lower=0)
    Ω_pd ∈ PDiagDomain(1)
    σ_add_pd ∈ RealDomain(lower=0)
  end

  @random begin
    ηpk ~ MvNormal(Ω_pk)
    ηpd ~ MvNormal(Ω_pd)
  end

  @pre begin
    CL = tvcl * exp(ηpk[1])
    Vc = tvvc * exp(ηpk[2])
    Q = tvq * exp(ηpk[3])
    Vp = tvvp * exp(ηpk[4])

    e_base = tvebase * exp(ηpd[1])
    ec_50 = tvec_50
    e_max = 1
    turn = tvturn
    k_out = 1/turn
    k_in0 = e_base*kout
  end

  @init begin
    Resp = e_base
  end

  @vars begin
    conc := Central/Vc
    e_drug := e_max*conc/(ec_50 + conc)
    k_in := k_in0*(1-e_drug)
  end

  @dynamics begin
   Central' = - (CL/Vc)*Central + (Q/Vp)*Peripheral - (Q/Vc)*Central
   Peripheral' = (Q/Vc)*Central - (Q/Vp)*Peripheral
   Resp' = k_in - k_out*Resp
  end

  @derived begin
    dv ~ @. Normal(conc, sqrt(conc^2*σ_prop_pk))
    resp ~ @. Normal(Resp, sqrt(σ_add_pd))
  end
end

If we want to extract LaTeX-formatted equations from the model then we first need to load Latexify.jl:

using Latexify

We can now latexify the different model blocks using calls like latexify(model, blockname) where blockname is a Symbol that indicates which model block to latexify. For example

latexify(model, :random)

\begin{align*} {\eta}pk &\sim \mathrm{MvNormal}\left( \Omega_{pk} \right) \\ {\eta}pd &\sim \mathrm{MvNormal}\left( \Omega_{pd} \right) \end{align*}

latexify(model, :pre) 

\begin{align*} CL &= tvcl \cdot e^{{\eta}pk_{1}} \\ Vc &= tvvc \cdot e^{{\eta}pk_{2}} \\ Q &= tvq \cdot e^{{\eta}pk_{3}} \\ Vp &= tvvp \cdot e^{{\eta}pk_{4}} \\ e_{base} &= tvebase \cdot e^{{\eta}pd_{1}} \\ ec_{50} &= tvec_{50} \\ e_{max} &= 1 \\ turn &= tvturn \\ k_{out} &= \frac{1}{turn} \\ k_{in0} &= e_{base} \cdot kout \end{align*}

latexify(model, :dynamics)

\begin{align*} \frac{dCentral(t)}{dt} =& \frac{Q \cdot Peripheral(t)}{Vp} - \frac{CL \cdot Central(t)}{Vc} - \frac{Q \cdot Central(t)}{Vc} \\ \frac{dPeripheral(t)}{dt} =& \frac{Q \cdot Central(t)}{Vc} - \frac{Q \cdot Peripheral(t)}{Vp} \\ \frac{dResp(t)}{dt} =& k_{in0} \cdot \left( 1 - \frac{e_{max} \cdot Central(t)}{\left( ec_{50} + \frac{Central(t)}{Vc} \right) \cdot Vc} \right) - k_{out} \cdot Resp(t) \end{align*}

latexify uses the information you input via the model macro, so some models contain information that others do not. Which blocknames are valid thus differ depending on the model that you use. In this example, the blocknames :param, :random, :pre, :init, :vars, :dynamics and :derived would be valid since these are the model blocks that we explicitly defined.

Using the latexify output

The latexify function returns a LaTeXString, which is pretty much a standard String with some overloads for making it display nicely. You can print this string to see the generated LaTeX code:

println(latexify(model, :dynamics))

\begin{align}
\frac{dCentral(t)}{dt} =& \frac{Q \cdot Peripheral(t)}{Vp} - \frac{CL \cdot Central(t)}{Vc} - \frac{Q \cdot Central(t)}{Vc} \\
\frac{dPeripheral(t)}{dt} =& \frac{Q \cdot Central(t)}{Vc} - \frac{Q \cdot Peripheral(t)}{Vp} \\
\frac{dResp(t)}{dt} =& k_{in0} \cdot \left( 1 - \frac{e_{max} \cdot Central(t)}{\left( ec_{50} + \frac{Central(t)}{Vc} \right) \cdot Vc} \right) - k_{out} \cdot Resp(t)
\end{align}

Different editors/environments have different support for more complex rendering. latexify automatically displays as nicely rendered equations in some environments, notably including notebooks. In other environments, like VSCode, you can call render to render the equations:

render(latexify(model, :dynamics))

alternatively

latexify(model, :dynamics) |> render

When we created Latexify.jl, we first thought that rendering would be somewhat frivolous, but quick access to readable equations has proven helpful during model development. Try for yourself.

Configurations

You can tune the latexify output to different keyword arguments. Many are provided directly from Latexify.jl, and some are specific to Pumas.jl models.

For example:

latexify(model, :dynamics; show_t=false, italicize=false, index=:bracket)

\begin{align*} \frac{\mathrm{dCentral}}{dt} =& \frac{Q \cdot \mathrm{Peripheral}}{Vp} - \frac{CL \cdot \mathrm{Central}}{Vc} - \frac{Q \cdot \mathrm{Central}}{Vc} \\ \frac{\mathrm{dPeripheral}}{dt} =& \frac{Q \cdot \mathrm{Central}}{Vc} - \frac{Q \cdot \mathrm{Peripheral}}{Vp} \\ \frac{\mathrm{dResp}}{dt} =& k_{in0} \cdot \left( 1 - \frac{e_{max} \cdot \mathrm{Central}}{\left( ec_{50} + \frac{\mathrm{Central}}{Vc} \right) \cdot Vc} \right) - k_{out} \cdot \mathrm{Resp} \end{align*}

For more options and features, have a look at Latexify.jl's documentation.

Stitching things together

Since the latexify output works well with other strings, you can quickly stitch different model parts together

model_representation = """
  \\section{Model dynamics}
  Here we might want to describe the model dynamics, as given by
  $(latexify(model, :dynamics)). 
  
  \\section{Random effects}
  The random effects are given by
  $(latexify(model, :random)).
  
  \\section{Conclusion}
  You have a lot of power to automate your model representation. 
  Don't forget to escape your backslashes. 
  """

println(model_representation)

\section{Model dynamics}
Here we might want to describe the model dynamics, as given by
\begin{align}
\frac{dCentral(t)}{dt} =& \frac{Q \cdot Peripheral(t)}{Vp} - \frac{CL \cdot Central(t)}{Vc} - \frac{Q \cdot Central(t)}{Vc} \\
\frac{dPeripheral(t)}{dt} =& \frac{Q \cdot Central(t)}{Vc} - \frac{Q \cdot Peripheral(t)}{Vp} \\
\frac{dResp(t)}{dt} =& k_{in0} \cdot \left( 1 - \frac{e_{max} \cdot Central(t)}{\left( ec_{50} + \frac{Central(t)}{Vc} \right) \cdot Vc} \right) - k_{out} \cdot Resp(t)
\end{align}
. 

\section{Random effects}
The random effects are given by
\begin{align}
{\eta}pk &\sim \mathrm{MvNormal}\left( \Omega_{pk} \right) \\
{\eta}pd &\sim \mathrm{MvNormal}\left( \Omega_{pd} \right)
\end{align}
.

\section{Conclusion}
You have a lot of power to automate your model representation. 
Don't forget to escape your backslashes. 

A note about future releases

The overall functionality described here is considered part of our public API and should thus not break during anything but a major release (Pumas.jl follow semantic versioning). However, only the intention behind a latexify command should be considered stable across non-breaking releases, not the exact formatting of the output string.