Analytical Model Solutions

Author

Patrick Kofod Mogensen

The topic of this tutorial is to introduce the use of analytical solutions in Pumas models. This includes practical considerations of how to use predefined models to obtain fast and accurate model solutions as part of simulation and estimation as well as a discussion of disadvantages and limitations.

1 Prerequisites

  1. Basic Julia Proficiency A fundamental understanding of the Julia programming language is beneficial.
  2. Pharmacometric Concepts A basic background in pharmacokinetics/pharmacodynamics (PK/PD) modeling or general NLME concepts ensures that terms such as ordinary differential equations, maximum likelihood estimation, and simulation are clear.
  3. Pumas Modeling Basic and Compartmental Modeling in Pumas The previous path and modules are assumed to have been completed such that the general workflow and the model components are known and understood.

2 Learning Goals

By the end of this tutorial, participants will be able to:

  • Efficiently use predefined compartmental models in Pumas
  • Understand the limitations of predefined models in practice
  • Write linear and affine ordinary differential equations to utilize the matrix exponential solver feature in Pumas

3 Analytical Solutions

3.1 Motivation and advantages

The goal of this tutorial is not to provide a theoretical background for analytical solutions to ordinary differential equations. The analytical solutions will be taken as given. However, it is beneficial to know why (and how) they are used.

3.1.1 Using Differential Equations for Drug Modeling

Researchers often use differential equations to understand how a drug moves through the body. The body is divided into separate “compartments,” such as the bloodstream or various tissues. The amount of drug in each compartment changes over time according to:

  1. Initial Conditions: At the start, there is a known initial amount of drug in each compartment.
  2. Transfer Rates: The drug may flow from one compartment to another, and it can also leave the body (for example, through the kidneys or liver). Each of these processes has a rate that describes how fast the transfer occurs.

Putting these together yields a system of differential equations that track how the drug amount changes in each compartment over time.

3.1.2 Linear, Affine, and Non-Linear Models

  • Linear Models: A model is called linear if the rate of change of the drug amount in a compartment depends in a strictly proportional way on the amounts in various compartments. If every compartment’s amount doubles, the rates of change also double.

  • Affine Models: An affine model is similar to a linear model but allows for a constant term. This constant term might represent a continuous infusion of the drug or some baseline production rate. In other words, an affine model can be viewed as a linear model plus some extra constant input or output.

  • Non-Linear Models: A model is considered non-linear if the rate of change does not scale in a simple proportional way with the drug amounts. Non-linearity can arise from processes such as enzyme saturation, cooperative binding, or complex feedback loops. In these cases, doubling the drug amount does not necessarily double the rate of change. Non-linear models often contain terms where variables multiply each other, appear in exponents, or otherwise go beyond simple proportional relationships.

3.1.3 Solving the Models: Eigenvalues and Matrix Exponentials

When a system of differential equations is linear (or affine), there are two main mathematical techniques that allow researchers to solve for the amounts of drug in each compartment in a fast and efficient way. Both rely on the fact that the solution to the system \frac{d\mathbf{x}(t)}{dt} = A \, \mathbf{x}(t), is given by \mathbf{x}(t) = e^{A t} \,\mathbf{x}(0), where A is the matrix encoding the rate constants and \mathbf{x}(0) represents the initial amounts of drug in each compartment. The question is then how to evaluate the matrix exponential function.

  1. Eigenproblem Approach: The matrix A represents the rates at which drug transfers between compartments. If we can diagonalize this matrix, it is much easier to evaluate the solution because the matrix exponential of a diagonal matrix is simply the usual exponential function of each diagonal element. By finding the eigenvalues and eigenvectors of this matrix, it is possible to build a closed-form solution for the drug amounts at any time. The eigenvalues give information about how fast different “modes” or “phases” of the system evolve.

  2. Matrix Exponential Approach: Another method involves evaluating the matrix exponential in an iterative way. This method is also fast and works in general for linear or affine systems. We will not go into details here, but it is essentially just an iterative method that is fast and generally accurate.

These two approaches are related and essentially yield the same solution. The former typically requires analysis up-front where as the latter can be used for any system we may provide using the differential equation interface in @dynamics.

3.1.4 Analytical vs. Numerical Solutions

  • Analytical (Closed-Form) Solutions: In many linear or affine cases, it is possible to find an explicit formula that describes exactly how the drug amounts evolve over time. Researchers can plug in any time value to get the drug amounts without resorting to iterative schemes. We will include the Matrix Exponential approach in this category because of its nice numerical properties.

  • Numerical Solutions: Non-linear models (and some more complicated cases) often do not have simple closed-form solutions. In these situations, we need to use more advanced solvers, typically using methods such as Euler’s method or Runge-Kutta. Numerical solutions can also be used when a closed-form solution would be too unwieldy to manipulate directly, even if one theoretically exists.

Further, there are models that are not linear or affine but still have analytical solutions such as Michaelis-Menten systems, but these are outside the scope of this tutorial.

3.2 Examples

Let us consider some examples.

using Pumas

3.2.1 A Linear ODE

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

The example above is a linear ordinary differential equation (linear ODE) because the right hand sides are sums of products between parameters and dynamical states. An analytical solution can be derived for this, and it can also be solved using the numerical matrix exponential method. Pumas can automatically detect if the matrix exponential is applicable by detecting if the right hand side is linear or affine.

linearode_model = @model begin
    @metadata begin
        desc = "Linear ODE Example"
        timeu = u"hr"
    end

    @param begin
        # PK parameters
        """
        Absorption Rate (L/h)
        """
        Ka ∈ RealDomain(lower = 0.0, init = 0.04)
        """
        Clearance (L/h)
        """
        CL ∈ RealDomain(lower = 0.0, init = 0.134)
        """
        Central Volume L
        """
        Vc ∈ RealDomain(lower = 0.0, init = 8.11)
    end

    @dynamics begin
        Depot' = -Ka * Depot
        Central' = Ka * Depot - CL / Vc * Central
    end
end
PumasModel
  Parameters: Ka, CL, Vc
  Random effects: 
  Covariates: 
  Dynamical system variables: Depot, Central
  Dynamical system type: Matrix exponential
  Derived: 
  Observed:

The output indicates Dynamical system type: Matrix exponential, showing that linearity was detected in this case and the fast solution method will be applied. To use the explicit analytical solution, the name of it in Pumas must be known. In this case, it is Depots1Central1 because there is one depot and one main compartment.

dept1cent1_model = @model begin
    @metadata begin
        desc = "Closed Form Oral Absorption Example"
        timeu = u"hr"
    end

    @param begin
        # PK parameters
        """
        Absorption Rate (L/h)
        """
        Ka ∈ RealDomain(lower = 0.0, init = 0.04)
        """
        Clearance (L/h)
        """
        CL ∈ RealDomain(lower = 0.0, init = 0.134)
        """
        Central Volume L
        """
        Vc ∈ RealDomain(lower = 0.0, init = 8.11)
    end

    @dynamics Depots1Central1
end
PumasModel
  Parameters: Ka, CL, Vc
  Random effects: 
  Covariates: 
  Dynamical system variables: Depot, Central
  Dynamical system type: Closed form
  Derived: 
  Observed:

The output indicates that a Closed form solution method is being used. Notice that closed form solutions require the parameters of the system (Ka, CL, Vc) to follow a specific naming scheme, and this scheme will always be available from the docstring ?Depots1Central1.

There are more models available in the documentation (external link). For example, there is a Depots2Central1 option for parallel absorption models but also more complicated models such as parent-metabolite models are available.

3.2.2 An Affine ODE

Another example would be an affine system.

@dynamics begin
    Central' = Kinfusion - CL / Vc * Central
end

Here Kinfusion represents a constant rate of infusion throughout the time course. If we define a minimal model

affineode_model = @model begin
    @metadata begin
        desc = "Affine ODE Example"
        timeu = u"hr"
    end

    @param begin
        # PK parameters
        """
        Infusion rate (L/h)
        """
        Kinfusion ∈ RealDomain(lower = 0.0, init = 0.04)
        """
        Clearance (L/h)
        """
        CL ∈ RealDomain(lower = 0.0, init = 0.134)
        """
        Central Volume L
        """
        Vc ∈ RealDomain(lower = 0.0, init = 8.11)
    end

    @dynamics begin
        Central' = Kinfusion - CL / Vc * Central
    end
end
PumasModel
  Parameters: Kinfusion, CL, Vc
  Random effects: 
  Covariates: 
  Dynamical system variables: Central
  Dynamical system type: Matrix exponential
  Derived: 
  Observed:

We once again see that the affine property of the system was detected and the solution method is indicated to be the matrix exponential.

Now, we typically hide this in the event system with an infusion event. If we re-write this to use the event system we can solve this model with a modified system where the infusion is implicitly handled through the event system.

@dynamics begin
    Central' = -CL / Vc * Central
end

Then we can also use the closed form model Central1.

infusion_closed_model = @model begin
    @metadata begin
        desc = "Affine ODE Example"
        timeu = u"hr"
    end

    @param begin
        """
        Clearance (L/h)
        """
        CL ∈ RealDomain(lower = 0.0, init = 0.134)
        """
        Central Volume L
        """
        Vc ∈ RealDomain(lower = 0.0, init = 8.11)
    end

    @dynamics Central1
end
PumasModel
  Parameters: CL, Vc
  Random effects: 
  Covariates: 
  Dynamical system variables: Central
  Dynamical system type: Closed form
  Derived: 
  Observed:

Here, it is clear that there will have to be a specification of events (in the analysis dataset for each individual) and dynamical evolution of the system conditional on the events (@dynamics). Even a constant rate that is not really an infusion can be modeled with an event and a model that then handles the distribution and elimination.

3.3 Limitations

Analytical (closed-form) solutions and the matrix exponential method are both fast and accurate, but they are not always applicable. Closed-form solutions can, in general, be derived for relatively large and complex models, although the expressions may become extremely cumbersome. In such situations, most of the efficiency can often be achieved by using the matrix exponential. However, two major limitations affect these two classes of model solutions:

  • They cannot include t (time) directly in any expressions because the methods rely on splitting the solution within each event interval.
  • They cannot handle partially closed-form and partially non-linear equations that are solved separately.

The first limitation can pose challenges for certain advanced models, for instance those involving sustained-release formulations with continuously changing rates over time. Nevertheless, in many practical cases, this is not a concern. Problems become more evident with the second limitation, which requires that all equations meet the assumptions. As an example, if there is oral absorption with a two-compartment distribution model (e.g., Depots1Central1Periph1) combined with an Emax PD model, closed-form solutions are not feasible because the PD model continuously depends on the PK model. Under these circumstances, the entire system must be solved numerically.

3.4 The Warfarin PKPD model

The Warfarin model that has been used in previous tutorials will be used again.

warfarin_pkpd_model = @model begin

    @metadata begin
        desc = "Warfarin PK/PD model"
        timeu = u"hr"
    end

    @param begin
        # PK parameters
        """
        Clearance (L/h/70kg)
        """
        pop_CL ∈ RealDomain(lower = 0.0, init = 0.134)
        """
        Central Volume L/70kg
        """
        pop_V ∈ RealDomain(lower = 0.0, init = 8.11)
        """
        Absorption time (h)
        """
        pop_tabs ∈ RealDomain(lower = 0.0, init = 0.523)
        """
        Lag time (h)
        """
        pop_lag ∈ RealDomain(lower = 0.0, init = 0.1)
        # PD parameters
        """
        Baseline
        """
        pop_e0 ∈ RealDomain(lower = 0.0, init = 100.0)
        """
        Emax
        """
        pop_emax ∈ RealDomain(init = -1.0)
        """
        EC50
        """
        pop_c50 ∈ RealDomain(lower = 0.0, init = 1.0)
        """
        Turnover
        """
        pop_tover ∈ RealDomain(lower = 0.0, init = 14.0)
        # Inter-individual variability
        """
          - ΩCL
          - ΩVc
          - ΩTabs
        """
        pk_Ω ∈ PDiagDomain([0.01, 0.01, 0.01])
        """
        Ωlag
        """
        lag_ω ∈ RealDomain(lower = 0.0, init = 0.1)
        """
          - Ωe0
          - Ωemax
          - Ωec50
          - Ωturn
        """
        pd_Ω ∈ PDiagDomain([0.01, 0.01, 0.01, 0.01])
        # Residual variability
        """
        Proportional residual error for drug concentration
        """
        σ_prop ∈ RealDomain(lower = 0.0, init = 0.00752)
        """
        Additive residual error for drug concentration (mg/L)
        """
        σ_add ∈ RealDomain(lower = 0.0, init = 0.0661)
        """
        Additive error for PCA
        """
        σ_fx ∈ RealDomain(lower = 0.0, init = 0.01)
    end

    @random begin
        # mean = 0, covariance = pk_Ω
        pk_η ~ MvNormal(pk_Ω)
        # mean = 0, standard deviation = lag_ω
        lag_η ~ Normal(0.0, lag_ω)
        # mean = 0, covariance = pd_Ω
        pd_η ~ MvNormal(pd_Ω)
    end

    @covariates begin
        """
        Scaling factor for Volume
        """
        FSZV
        """
        Scaling factor for Clearance
        """
        FSZCL
    end

    @pre begin
        # PK
        CL = FSZCL * pop_CL * exp(pk_η[1])
        Vc = FSZV * pop_V * exp(pk_η[2])
        tabs = pop_tabs * exp(pk_η[3])
        Ka = log(2) / tabs
        # PD
        e0 = pop_e0 * exp(pd_η[1])
        emax = pop_emax * exp(pd_η[2])
        c50 = pop_c50 * exp(pd_η[3])
        tover = pop_tover * exp(pd_η[4])
        kout = log(2) / tover
        rin = e0 * kout
    end

    @dosecontrol begin
        lags = (Depot = pop_lag * exp(lag_η),)
    end

    @init begin
        Turnover = e0
    end

    @vars begin
        cp := Central / Vc
        ratein := Ka * Depot
        pd := 1 + emax * cp / (c50 + cp)
    end

    @dynamics begin
        Depot' = -ratein
        Central' = ratein - CL * cp
        Turnover' = rin * pd - kout * Turnover
    end

    @derived begin
        """
        Warfarin Concentration (mg/L)
        """
        conc ~ @. Normal(cp, sqrt((σ_prop * cp)^2 + σ_add^2))
        """
        PCA
        """
        pca ~ @. Normal(Turnover, σ_fx)
    end
end
PumasModel
  Parameters: pop_CL, pop_V, pop_tabs, pop_lag, pop_e0, pop_emax, pop_c50, pop_tover, pk_Ω, lag_ω, pd_Ω, σ_prop, σ_add, σ_fx
  Random effects: pk_η, lag_η, pd_η
  Covariates: FSZV, FSZCL
  Dynamical system variables: Depot, Central, Turnover
  Dynamical system type: Nonlinear ODE
  Derived: conc, pca
  Observed: conc, pca

By the looks of it, this is a linear model but looks can be deceiving!

@dynamics begin
    Depot' = -ratein
    Central' = ratein - CL * cp
    Turnover' = rin * pd - kout * Turnover
end

Everything on the right hand sides are nothing but additions (subtractions) and multiplications, but one has to remember that aliases were created in the @vars block. The problem becomes clear on substituting the expressions.

@dynamics begin
    Depot' = -Ka * Depot
    Central' = Ka * Depot - CL * Central / Vc
    Turnover' = rin * (1 + emax * Central / Vc / (c50 + Central / Vc)) - kout * Turnover
end

The problem is that Central appears in the denominator of the rin part of the Turnover model. Looking at the model definition output, the printed information says Dynamical system type: Nonlinear ODE that confirms the suspicion.

It is not the case that because models are non-linear they cannot have an analytical solution, but Pumas does not currently have a closed form model for the above model.

4 Summary

In summary, this discussion focuses on how to leverage predefined analytical solutions and a matrix exponential solver in Pumas. The main advantages and disadvantages of using analytical solutions were reviewed, and an assessment was conducted to determine whether these methods could be applied to the Warfarin model used extensively in this series. It was not possible to use analytical solutions, even though the absorption, distribution, and elimination components fit the Depots1Central1 assumptions, because the PD model includes a term with Central in the denominator, making the overall model non-linear.