Dynamical Systems

using Pumas
using Pumas.Catalyst
using Pumas.ModelingToolkit
using AlgebraOfGraphics
using CairoMakie

using Random

In this tutorial, we will discuss different approaches for defining dynamical systems in a Pumas model. These dynamical systems are differential equation models that describe the (typically non-linear) function from the parameters to the derived variables. They are the defined in the @dynamics or alternatively the @reactions block of a Pumas model.

We compare the different approaches by implementing a simple one-compartment model in 9 different ways.

1 Closed-Form Models (Pumas)

Some compartment models with analytical solutions are predefined in Pumas. They are referred to by name and their parameters have to be defined in the @pre block. Since the solutions of these models are available in closed form, simulation and estimation can be much faster than for generic differential equation model.

For instance, Depots1Central1 is a one-compartment model with a central compartment Central and a depot Depot. More formally, the model is defined as

\[ \begin{aligned} \frac{\mathrm{d}}{\mathrm{d}t} \operatorname{Depot}(t) &= - \operatorname{Ka} \operatorname{Depot}(t), \\ \frac{\mathrm{d}}{\mathrm{d}t} \operatorname{Central}(t) &= \operatorname{Ka} \operatorname{Depot}(t) - \frac{\operatorname{CL}}{\operatorname{Vc}} \operatorname{Central}(t), \end{aligned} \]

with positive parameters \(\operatorname{Ka}\) (absorption rate), \(\operatorname{CL}\) (clearance), and \(\operatorname{Vc}\) (volume of distribution). Its dynamics can be used in a Pumas model by adding

@dynamics Depots1Central1

inside of a @model or @emmodel definition.

We define a Pumas model whose dynamics are governed by the one-compartment model Depots1Central1. For illustration purposes, we omit fixed and random effects, covariates and dependent variables. In addition to the @dynamics block only a @pre block with the definition of the ODE parameters Ka, CL and Vc is required. To simplify the comparison with the other models defined below, we also define a non-zero initial value of Depot.

closedform = @model begin
    @init begin
        Depot = 2.5
    end
    @pre begin
        Ka = 3.8
        CL = 3.2
        Vc = 16.4
    end
    @dynamics Depots1Central1
end

PumasModel
  Parameters:
  Random effects:
  Covariates:
  Dynamical system variables: Depot, Central
  Dynamical system type: Closed form
  Derived:
  Observed:

The model summary shows that indeed this Pumas model uses a “Closed form” dynamical system with state variables Depot and Central.

2 Differential Equations (Pumas)

Generic differential equations can be defined in the @dynamics block by directly stating the differential equations as a system of equations of the form

x' = ...

with a derivative (specified by the suffix ') of a variable x on the left-hand side and the symbolic expression of the derivative on the right-hand side.

Tip

You can use time variable t and variables defined in the @vars block on the right-hand side of the derivative expressions.

For instance, the one-compartment model Depots1Central1 can be declared as

@dynamics begin
    Depot' = -Ka * Depot
    Central' = Ka * Depot - (CL / Vc) * Central
end

By default, Pumas checks whether the derivative expressions define a linear differential equation system, i.e., a differential equation of the form

\[ \frac{\mathrm{d}}{\mathrm{d}t} \begin{bmatrix} x_1(t) \\ x_2(t) \\ \vdots \\ x_n(t) \\ \end{bmatrix} = A \begin{bmatrix} x_1(t) \\ x_2(t) \\ \vdots \\ x_n(t) \\ \end{bmatrix} + b \]

where matrix \(A \in \mathbb{R}^{n \times n}\) and vector \(b \in \mathbb{R}^n\) are constant. If a linear differential equation is detected, then Pumas simulates it by evaluating the formal solution

\[ \begin{bmatrix} x_1(t) \\ x_2(t) \\ \vdots \\ x_n(t) \\ \end{bmatrix} = \big(\exp{(At)} - I_n \big) A^{-1} b + \exp{(At)} \begin{bmatrix} x_1(0) \\ x_2(0) \\ \vdots \\ x_n(0) \\ \end{bmatrix}, \]

assuming that \(A\) is invertible.

Linearity detection can be disabled with the checklinear = false setting in the @options block.

With Linearity Detection

linear = @model begin
    @init begin
        Depot = 2.5
    end
    @pre begin
        Ka = 3.8
        CL = 3.2
        Vc = 16.4
    end
    @dynamics begin
        Depot' = -Ka * Depot
        Central' = Ka * Depot - (CL / Vc) * Central
    end
end

PumasModel
  Parameters:
  Random effects:
  Covariates:
  Dynamical system variables: Depot, Central
  Dynamical system type: Matrix exponential
  Derived:
  Observed:

The description “Dynamical system type: Matrix exponential” in the summary of the Pumas model indicates that Pumas correctly discovered that the differential equations define a linear system with an explicit solution.

Without Linearity Detection

generic = @model begin
    @options begin
        checklinear = false
    end
    @init begin
        Depot = 2.5
    end
    @pre begin
        Ka = 3.8
        CL = 3.2
        Vc = 16.4
    end
    @dynamics begin
        Depot' = -Ka * Depot
        Central' = Ka * Depot - (CL / Vc) * Central
    end
end

PumasModel
  Parameters:
  Random effects:
  Covariates:
  Dynamical system variables: Depot, Central
  Dynamical system type: Nonlinear ODE
  Derived:
  Observed:

The output “Dynamical system type: Nonlinear ODE” shows that the dynamics of this Pumas model are defined as a generic differential equation problem that is simulated without exploiting linearity.

3 Differential Equations (ModelingToolkit)

The definition of the dynamical system can be decoupled from the definition of the Pumas model. A symbolic differential equation system sys defined with ModelingToolkit.jl can be used as dynamical system of a Pumas model by specifying

@dynamics sys

inside of the @model or @emmodel macro.

This separation of the Pumas model and its dynamical system can be useful if the dynamical system is of interest in itself and one wants to study some of its properties such as stability, bifurcations, steady states, parameter sensitivity or structural identifiability. It also makes it easy to re-use the same dynamics in different models, similar to closed-form dynamics such as Depots1Central1. Additionally, it allows you to import dynamical systems from external sources, such as SBML or CellML files.

We define a symbolic ModelingToolkit.ODESystem that corresponds to the one-compartment model Depots1Central1. For more information about how to define and work with such symbolic systems, please check the documentation of ModelingToolkit.jl.

# Independent variable
@independent_variables t

# Differential operator
D = Differential(t)

# Definition of the dynamical system
@mtkmodel MTKDepots1Central1 begin
    # Parameters
    @parameters begin
        Ka = 3.8 # default value
        CL
        Vc
    end

    # State variables (with and)
    @variables begin
        Depot(t) = 2.5 # initial value
        Central(t)
    end

    # Differential equations
    @equations begin
        D(Depot) ~ -Ka * Depot
        D(Central) ~ Ka * Depot - (CL / Vc) * Central
    end
end

# Instantiated system
@mtkbuild mtksys = MTKDepots1Central1()

Model mtksys:
Equations (2):
  2 standard: see equations(mtksys)
Unknowns (2): see unknowns(mtksys)
  Central(t)
  Depot(t) [defaults to 2.5]
Parameters (3): see parameters(mtksys)
  Ka [defaults to 3.8]
  Vc
  CL

We use this dynamical system to re-define the one-compartment model. This time, however, we do not define the absorption rate parameter Ka in the @pre block. This implies that in the Pumas model Ka is fixed to its default value of 3.8 specified in the symbolic system.

With Linearity Detection

mtk_linear = @model begin
    @pre begin
        CL = 3.2
        Vc = 16.4
    end
    @dynamics mtksys
end

┌ Warning: The model uses the following non-zero initial values defined implicitly in the dynamical system:
│   Depot = 2.5
└ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/6G31F/src/dsl/model_macro.jl:1442

PumasModel
  Parameters:
  Random effects:
  Covariates:
  Dynamical system variables: Central, Depot
  Dynamical system type: Matrix exponential
  Derived:
  Observed:

The model summary shows that Pumas correctly detected that the dynamical system defines a linear differential equation (“Matrix exponential”).

Without Linearity Detection

Again we can disable the linearity check with the checklinear = false option.

mtk_generic = @model begin
    @options begin
        checklinear = false
    end
    @init begin
        Depot = 2.5
    end
    @pre begin
        CL = 3.2
        Vc = 16.4
    end
    @dynamics mtksys
end

PumasModel
  Parameters:
  Random effects:
  Covariates:
  Dynamical system variables: Central, Depot
  Dynamical system type: Nonlinear ODE
  Derived:
  Observed:

Pumas Warnings

When defining mtk_linear, Pumas warns us that the initial value of Depot is 2.5, as specified in the symbolic system mtksys:

ModelingToolkit.defaults(mtksys)

Dict{Any, Any} with 6 entries:
  Initial(Centralˍt(t)) => false
  Ka                    => 3.8
  Initial(Depotˍt(t))   => false
  Initial(Depot(t))     => false
  Depot(t)              => 2.5
  Initial(Central(t))   => false

However, Pumas doesn’t show a warning when defining mtk_generic since an initial value of Depot is defined in its @init block, and hence the Pumas model does not use any initial values of the symbolic system.

4 Reactions (Pumas)

Instead of specifying the differential equations of the dynamical system in a @dynamics block inside of the @model or @emmodel macro, we can specify the transitions (“reactions”) between the state variables/compartments in a @reactions block. Arguably, stating the transitions can be both more convenient and clearer, in particular with increasing number of variables.

The transitions have to be defined with the Catalyst.jl syntax. In each line of the @reactions block a transition is defined. The most basic transition is

k, A --> B

which describes a transition of rate k from state A to state B.

Our one-compartment model consists of two transitions:

• The drug transitions from Depot to Central with absorption rate Ka

• The drug is eliminated from Central with rate CL / Vc

In the Catalyst.jl syntax, elimination from the system is described as a transition to a 0 state. Put together, the one-compartment model can be reimplemented as:

With Linearity Detection

reactions_linear = @model begin
    @init begin
        Depot = 2.5
    end
    @pre begin
        Ka = 3.8
        CL = 3.2
        Vc = 16.4
    end
    @reactions begin
        Ka, Depot --> Central
        CL / Vc, Central --> 0
    end
end

PumasModel
  Parameters:
  Random effects:
  Covariates:
  Dynamical system variables: Depot, Central
  Dynamical system type: Matrix exponential
  Derived:
  Observed:

The model summary shows that also in this case Pumas correctly detects that the system describes a linear differential equation.

Without Linearity Detection

reactions_generic = @model begin
    @options begin
        checklinear = false
    end
    @init begin
        Depot = 2.5
    end
    @pre begin
        Ka = 3.8
        CL = 3.2
        Vc = 16.4
    end
    @reactions begin
        Ka, Depot --> Central
        CL / Vc, Central --> 0
    end
end

PumasModel
  Parameters:
  Random effects:
  Covariates:
  Dynamical system variables: Depot, Central
  Dynamical system type: Nonlinear ODE
  Derived:
  Observed:

Exclusive Blocks

It is not allowed for a model to contain both a @dynamics and an @reactions block.

@model begin
    @dynamics Depots1Central1
    @reactions begin
        Ka2, Depot --> Central2
    end
end

LoadError: ArgumentError: Both `@reactions` and `@dynamics` blocks cannot be defined in the same model.
in expression starting at /build/run/_work/PumasTutorials.jl/PumasTutorials.jl/tutorials/introduction/dynamical_systems.qmd:393
Stacktrace:
  [1] check_double_def!(blocklist::Set{Symbol}, x::Symbol)
    @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/6G31F/src/dsl/model_macro.jl:2600
  [2] _model(expr::Expr, currmodule::Module, source::LineNumberNode)
    @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/6G31F/src/dsl/model_macro.jl:3044
  [3] var"@model"(__source__::LineNumberNode, __module__::Module, ex::Any)
    @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/6G31F/src/dsl/model_macro.jl:3425
  [4] eval
    @ ./boot.jl:430 [inlined]
  [5] (::QuartoNotebookWorker.var"#15#16"{Module, Expr})()
    @ QuartoNotebookWorker ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/QuartoNotebookRunner/fVEXT/src/QuartoNotebookWorker/src/render.jl:222
  [6] (::QuartoNotebookWorker.Packages.IOCapture.var"#5#9"{DataType, QuartoNotebookWorker.var"#15#16"{Module, Expr}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}})()
    @ QuartoNotebookWorker.Packages.IOCapture ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/QuartoNotebookRunner/fVEXT/src/QuartoNotebookWorker/src/vendor/IOCapture/src/IOCapture.jl:170
  [7] with_logstate(f::QuartoNotebookWorker.Packages.IOCapture.var"#5#9"{DataType, QuartoNotebookWorker.var"#15#16"{Module, Expr}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}}, logstate::Base.CoreLogging.LogState)
    @ Base.CoreLogging ./logging/logging.jl:524
  [8] with_logger(f::Function, logger::Base.CoreLogging.ConsoleLogger)
    @ Base.CoreLogging ./logging/logging.jl:635
  [9] capture(f::QuartoNotebookWorker.var"#15#16"{Module, Expr}; rethrow::Type, color::Bool, passthrough::Bool, capture_buffer::IOBuffer, io_context::Vector{Pair{Symbol, Any}})
    @ QuartoNotebookWorker.Packages.IOCapture ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/QuartoNotebookRunner/fVEXT/src/QuartoNotebookWorker/src/vendor/IOCapture/src/IOCapture.jl:167
 [10] capture
    @ ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/QuartoNotebookRunner/fVEXT/src/QuartoNotebookWorker/src/render.jl:248 [inlined]
 [11] io_capture(f::Function; cell_options::Dict{Any, Any}, kws::@Kwargs{rethrow::DataType, color::Bool, io_context::Vector{Pair{Symbol, Any}}})
    @ QuartoNotebookWorker ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/QuartoNotebookRunner/fVEXT/src/QuartoNotebookWorker/src/render.jl:250
 [12] io_capture
    @ ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/QuartoNotebookRunner/fVEXT/src/QuartoNotebookWorker/src/render.jl:246 [inlined]
 [13] include_str(mod::Module, code::String; file::String, line::Int64, cell_options::Dict{Any, Any})
    @ QuartoNotebookWorker ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/QuartoNotebookRunner/fVEXT/src/QuartoNotebookWorker/src/render.jl:201
 [14] #invokelatest#2
    @ ./essentials.jl:1057 [inlined]
 [15] invokelatest
    @ ./essentials.jl:1052 [inlined]
 [16] #6
    @ ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/QuartoNotebookRunner/fVEXT/src/QuartoNotebookWorker/src/render.jl:18 [inlined]
 [17] with_inline_display(f::QuartoNotebookWorker.var"#6#7"{String, String, Int64, Dict{Any, Any}}, cell_options::Dict{Any, Any})
    @ QuartoNotebookWorker ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/QuartoNotebookRunner/fVEXT/src/QuartoNotebookWorker/src/InlineDisplay.jl:31
 [18] _render_thunk(thunk::Function, code::String, cell_options::Dict{Any, Any}, is_expansion_ref::Base.RefValue{Bool}; inline::Bool)
    @ QuartoNotebookWorker ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/QuartoNotebookRunner/fVEXT/src/QuartoNotebookWorker/src/render.jl:43
 [19] _render_thunk
    @ ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/QuartoNotebookRunner/fVEXT/src/QuartoNotebookWorker/src/render.jl:35 [inlined]
 [20] #render#5
    @ ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/QuartoNotebookRunner/fVEXT/src/QuartoNotebookWorker/src/render.jl:15 [inlined]
 [21] render(code::String, file::String, line::Int64, cell_options::Dict{Any, Any})
    @ QuartoNotebookWorker ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/QuartoNotebookRunner/fVEXT/src/QuartoNotebookWorker/src/render.jl:1
 [22] render(::String, ::Vararg{Any}; kwargs::@Kwargs{})
    @ Main ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/QuartoNotebookRunner/fVEXT/src/startup.jl:145
 [23] top-level scope
    @ none:1

5 Reactions (Catalyst)

Reactions can be decoupled from the Pumas model as well. In fact, a Catalyst.ReactionSystem is another example of a symbolic ModelingToolkit system.

Since the @reactions blocks already use Catalyst.jl syntax, arguably defining a separate Catalyst.ReactionSystem is even simpler than defining a ModelingToolkit.ODESystem such as the MTKDepots1Central1 example above. Basically, we only have to replace @reactions with Catalyst.@reaction_networks.

For instance, the one-compartment model Depots1Central corresponds to the following reaction network:

@reaction_network CatalystDepots1Central1 begin
    Ka, Depot --> Central
    CL / Vc, Central --> 0
end

Model CatalystDepots1Central1:
Unknowns (2): see unknowns(CatalystDepots1Central1)
  Depot(t)
  Central(t)
Parameters (3): see parameters(CatalystDepots1Central1)
  Ka
  CL
  Vc

If desired, we can define default parameter values and initial values of the states by adding @parameters and @species sections inside of the @reaction_network:

catalystsys = @reaction_network CatalystDepots1Central1 begin
    @parameters Ka = 3.8
    @species Depot(t) = 2
    Ka, Depot --> Central
    CL / Vc, Central --> 0
end
catalystsys = complete(catalystsys) # only required in Catalyst < 0.14

Model CatalystDepots1Central1:
Unknowns (2): see unknowns(CatalystDepots1Central1)
  Depot(t) [defaults to 2]
  Central(t)
Parameters (3): see parameters(CatalystDepots1Central1)
  Ka [defaults to 3.8]
  CL
  Vc

With Linearity Detection

catalyst_linear = @model begin
    @init begin
        Depot = 2.5
    end
    @pre begin
        CL = 3.2
        Vc = 16.4
    end
    @dynamics catalystsys
end

PumasModel
  Parameters:
  Random effects:
  Covariates:
  Dynamical system variables: Depot, Central
  Dynamical system type: Matrix exponential
  Derived:
  Observed:

The model summary confirms that Pumas detected successfully that the reaction network defines a linear differential equation.

Without Linearity Detection

catalyst_generic = @model begin
    @options begin
        checklinear = false
    end
    @init begin
        Depot = 2.5
    end
    @pre begin
        CL = 3.2
        Vc = 16.4
    end
    @dynamics catalystsys
end

PumasModel
  Parameters:
  Random effects:
  Covariates:
  Dynamical system variables: Depot, Central
  Dynamical system type: Nonlinear ODE
  Derived:
  Observed:

As expected, we see that Pumas treats the system as a generic differential equation after disabling the linearity check.

6 Conclusion

In this tutorial, we discussed multiple approaches for defining the dynamical system of a Pumas model. Figure 1 illustrates that all models presented above yield the same simulations.

let
    # Perform simulations
    models = (
        ("closedform", closedform),
        ("linear", linear),
        ("generic", generic),
        ("mtk_linear", mtk_linear),
        ("mtk_generic", mtk_generic),
        ("reactions_linear", reactions_linear),
        ("reactions_generic", reactions_generic),
        ("catalyst_linear", catalyst_linear),
        ("catalyst_generic", catalyst_generic),
    )
    sims = mapreduce(vcat, models) do (name, model)
        sim = simobs(model, Subject(), (;); obstimes = range(0, 10; length = 1001))
        return DataFrame(; sim.time, sim.dynamics.Central, name)
    end

    # Plot simulations of Central
    plt = data(sims) * mapping(:time, :Central, layout = :name) * visual(Lines)
    draw(plt)
end

Figure 1: Simulations from the models in this tutorial.

If a closed-form model of the dynamical system is available in Pumas, probably you should use it: Closed-form models tend to be fast both in simulation and estimation. If you want to study a dynamical system in more detail, share it between Pumas models, or import it from a file, then you should consider defining it separately from the Pumas model.