Distributed Delay Absorption

The traditional way to model absorption delay is to use a lag time parameter (e.g., ALAG parameter in NONMEM). The lag time shifts the time of dosing as if the drug was in fact administered at a later time. Introducing a drug into the system after a certain lag time indicates a sudden rise in the absorption rate from an initial value of zero. However, this approach is quite non-physiological. Furthermore, the discontinuous nature of the resulting concentration-time curve can pose challenges when attempting to determine the optimal parameter values of the model through search algorithms. This can happen because the objective function can be non-differentiable at the lag time, and arises on situations where the lag or duration parameter is modeled with between-subject variability. One suggestion to circumvent this issue is to avoid modeling approaches that introduce non-differentiabilities in the likelihood, e.g., by using delay differential equation methods such as transit compartment models [14].

The use of transit compartment models to describe delays in onset of absorption was developed by [15], and the main idea is to estimate from data the optimal number of transit compartments and the rate of mass transfer between these compartments. This is achieved by using a continuous distribution for the number of transit compartments instead of using a discrete number of transit compartments. Hence, the number of transit compartments is not fixed, but rather a random variable.

In Pumas, the distributed delay absorption model is implemented using the @delay macro. The macro @delay takes two arguments: a distribution and a compartment. distribution can be any of the distributions supported by Pumas listed in Defining NLME models in Pumas. compartment must be one of the compartments defined in the @dynamics block, e.g. Central. It returns a real number that represents the estimated number of transit compartments. This macro should be used in the @pre block to compute an absorption rate explicitly added to the ODE in the @dynamics block.

Below are two examples: a Weibull absorption model of [16]; and a Gamma absorption model of [17].

@delay(
    distribution::Distribution,
    cmt,
) -> Real

mdl_weibull = @model begin
  @param begin
    θCb ∈ RealDomain(lower = 0.0)
    θVc ∈ RealDomain(lower = 0.0)
    θVm ∈ RealDomain(lower = 10.0)
    θKm ∈ RealDomain(lower = 0.0)

    θA ∈ RealDomain(lower = 0.0)
    θB ∈ RealDomain(lower = 0.0)

    θf ∈ RealDomain(lower = 0.0)

    σ ∈ RealDomain(lower = 0.0)
  end

  @pre begin
    Baseline = θCb * θVc
    Vc = θVc
    VVm = Vc * θVm
    VKm = Vc * θKm

    λ = θA^(-1 / θB)
    B = θB

    f = θf

    dose_density = @delay(Weibull(B, λ), Central)
  end

  @dosecontrol begin
    bioav = (; Central = 0.0)
  end

  @init begin
    Central = Baseline
  end

  @dynamics begin
    Central' = VVm * (Baseline / (VKm + Baseline) - Central / (VKm + Central)) + f * Vc * dose_density
  end

  @derived begin
    μ := @. Central ./ Vc
    y ~ @. Gamma(1 / σ^2, μ * σ^2)
  end
end

mdl_gamma_delay = @model begin
  @param begin
    θn ∈ RealDomain(lower = 0.0)
    θmat ∈ RealDomain(lower = 0.0)
    θCL ∈ RealDomain(lower = 0.0)
    θVc ∈ RealDomain(lower = 0.0)

    σ ∈ RealDomain(lower = 0.0)
  end

  @pre begin
    kₑ = θCL / θVc
    Vc = θVc
    dose_density = @delay(Gamma(θn, θmat / θn), Central)
  end

  @dosecontrol begin
    bioav = (; Central = 0.0)
  end

  @dynamics begin
    Central' = dose_density - kₑ * Central
  end

  @derived begin
    μ := Central ./ Vc
    y ~ @. Gamma(1 / σ^2, μ * σ^2)
  end
end