Absorption Models

Authors

Vijay Ivaturi

Jose Storopoli

Patrick Kofod Mogensen

using Random
using Pumas
using PumasUtilities
using CairoMakie

Random.seed!(1234)  # Set random seed for reproducibility

1 Introduction

Absorption modeling is an integral part of pharmacokinetic modeling for drugs delivered by any route of administration other than intravenous (IV) (e.g. oral, subcutaneous, intramuscular, transdermal, nasal, etc.). In the simplest case, drugs administered orally may undergo first-order absorption, whereby the absorption from the depot compartment (i.e. gut) to the central compartment (i.e. plasma) is a first-order process. However, many more complex situations arise in practice, some of which will becovered in this tutorial.

We will cover the following absorption models:

First-order
Zero-order
Parallel zero-order and first-order
Two parallel first-order processes
Transit compartments and absorption according to the Erlang distribution
Absorption according to the Gamma distribution
Absorption according to the Weibull distribution

This tutorial is meant to give a broad overview of how to handle different absorption models in Pumas and as a reference. Other tutorials will go into details about specific models and case-studies. We restrict ourselves to a single subject without any random effects (by omitting the @random block) or observational noise (by setting σ=0.0).

For a more in-depth discussion of simulating data with Pumas, please see the tutorial generating and simulation populations.

2 First-Order Absorption

Note

Note that we have written the model using differential equations, but we could equally have specified an “analytical” (closed-form) solution using this @dynamics block instead (see documentation on analytical solutions):

@dynamics Depots1Central1

Here is the Pumas code for a 1-compartment model with first-order absorption, diagonal random-effects structure, and allometric scaling on clearance and volume:

foabs = @model begin

    @param begin
        tvcl ∈ RealDomain(; lower = 0)
        tvvc ∈ RealDomain(; lower = 0)
        tvka ∈ RealDomain(; lower = 0)
        σ ∈ RealDomain(; lower = 0)
    end

    @pre begin
        CL = tvcl
        Vc = tvvc
        Ka = tvka
    end

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

    @derived begin
        conc = @. Normal(Central / Vc, σ)
    end
end

PumasModel
  Parameters: tvcl, tvvc, tvka, σ
  Random effects: 
  Covariates: 
  Dynamical system variables: Depot, Central
  Dynamical system type: Matrix exponential
  Derived: conc
  Observed: conc

param1 = (; tvcl = 5, tvvc = 20, tvka = 1, σ = 0.0)

dose1 = DosageRegimen(100; time = 0)
subj = Subject(; id = "1", events = dose1, observations = (; conc = nothing))

sims1 = simobs(foabs, subj, param1; obstimes = 0:0.1:48)

SimulatedObservations
  Simulated variables: conc
  Time: 0.0:0.1:48.0

sim_plot(
    sims1;
    observations = [:conc],
    figure = (; fontsize = 18),
    axis = (;
        xlabel = "Time (hr)",
        ylabel = "Predicted Concentration (ng/mL)",
        title = "First-order absorption",
    ),
)

3 Zero-Order Absorption

Note

Notice that there is no depot compartment, and that instead the absorption takes place in the @dosecontrol block by specifying a duration for the zero-order process. In order for this to work, we need to set the rate data item in DosageRegimen to the value -2; this is a clue to the Pumas engine that the rate should be derived from the duration and amt.

Zero-order absorption is less common than first-order absorption. It is essentially like an IV infusion, but where the duration of the infusion is an estimated parameter rather than a known quantity. Here is the Pumas code for the same 1-compartment model but with zero-order absorption:

zoabs = @model begin
    @param begin
        tvcl ∈ RealDomain(; lower = 0)
        tvvc ∈ RealDomain(; lower = 0)
        tvdur ∈ RealDomain(; lower = 0)
        σ ∈ RealDomain(; lower = 0)
    end

    @pre begin
        CL = tvcl
        Vc = tvvc
    end

    @dosecontrol begin
        duration = (; Central = tvdur)
    end

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

    @derived begin
        conc = @. Normal(Central / Vc, σ)
    end

end

PumasModel
  Parameters: tvcl, tvvc, tvdur, σ
  Random effects: 
  Covariates: 
  Dynamical system variables: Central
  Dynamical system type: Matrix exponential
  Derived: conc
  Observed: conc

param2 = (; tvcl = 0.792, tvvc = 13.7, tvdur = 5.0, σ = 0.0)
dose2 = DosageRegimen(
    100;
    time = 0,
    rate = -2, # Note: rate must be -2 for duration modeling
)
subj2 = map(i -> Subject(; id = i, events = dose2, observations = (; conc = nothing)), 1:10)
sims2 = simobs(zoabs, subj2, param2; obstimes = 0:0.1:48)

Simulated population (Vector{<:Subject})
  Simulated subjects: 10
  Simulated variables: conc

sim_plot(
    zoabs,
    sims2;
    observations = [:conc],
    figure = (; fontsize = 18),
    axis = (;
        xlabel = "Time (hr)",
        ylabel = "Predicted Concentration (ng/mL)",
        title = "Zero-order absorption",
    ),
)

4 Parallel Zero-Order and First-Order Absorption

This is a more complex absorption model, with both a zero-order process and a first-order process operating simultaneously. Essentially, the total dose is split into two parts, one of which undergoes first-order absorption, and the other zero-order absorption (how this splitting is typically accomplished from a technical perspective is to treat each actual dose as two virtual doses with different bioavailable fractions, as in the example below). Furthermore, we will assume that the first-order process is subject to a lag (i.e., the absorption begins after a certain time delay, which is a model parameter).

zofo_paral_abs = @model begin

    @param begin
        tvcl ∈ RealDomain(; lower = 0)
        tvvc ∈ RealDomain(; lower = 0)
        tvka ∈ RealDomain(; lower = 0)
        tvdur ∈ RealDomain(; lower = 0)
        tvbio ∈ RealDomain(; lower = 0)
        tvlag ∈ RealDomain(; lower = 0)
        σ ∈ RealDomain(; lower = 0)
    end

    @pre begin
        CL = tvcl
        Vc = tvvc
        Ka = tvka
    end

    @dosecontrol begin
        duration = (; Central = tvdur)
        bioav = (; Depot = tvbio, Central = 1 - tvbio)
        lags = (; Depot = tvlag)
    end

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

    @derived begin
        conc = @. Normal(Central / Vc, σ)
    end

end

PumasModel
  Parameters: tvcl, tvvc, tvka, tvdur, tvbio, tvlag, σ
  Random effects: 
  Covariates: 
  Dynamical system variables: Depot, Central
  Dynamical system type: Matrix exponential
  Derived: conc
  Observed: conc

Note

With these parameter values, the first-order process starts 1 hour after administration, while the zero-order process lasts 2 hours. Hence, there is one hour where they overlap.

param3 = (; tvcl = 5, tvvc = 50, tvka = 1.2, tvdur = 2, tvbio = 0.5, tvlag = 1, σ = 0.0)
dose_zo = DosageRegimen(
    100;
    time = 0,
    cmt = 2,
    rate = -2, # Note: rate must be -2 for duration modeling
    evid = 1,
)
dose_fo = DosageRegimen(100; time = 0, cmt = 1, rate = 0, evid = 1)
dose3 = DosageRegimen(dose_fo, dose_zo)  # Actual dose is made up of 2 virtual doses
subj3 = map(i -> Subject(; id = i, events = dose3, observations = (; conc = nothing)), 1:10)
sims3 = simobs(zofo_paral_abs, subj3, param3; obstimes = 0:0.1:48)

Simulated population (Vector{<:Subject})
  Simulated subjects: 10
  Simulated variables: conc

sim_plot(
    sims3;
    observations = [:conc],
    figure = (; fontsize = 18),
    axis = (;
        xlabel = "Time (hr)",
        ylabel = "Predicted Concentration (ng/mL)",
        title = "Zero- and first-order parallel absorption",
    ),
)

5 Absorption by Two Parallel First-Order Processes

This is similar to the example above, except that both parallel processes are first-order rather than one being zero-order. We will refer to the two processes as “immediate release” (IR) and “slow release” (SR); only the SR process is subject to a lag. The Pumas code is as follows:

two_parallel_foabs = @model begin

    @param begin
        tvcl ∈ RealDomain(; lower = 0)
        tvvc ∈ RealDomain(; lower = 0)
        tvka1 ∈ RealDomain(; lower = 0)
        tvka2 ∈ RealDomain(; lower = 0)
        tvlag ∈ RealDomain(; lower = 0)
        tvbio ∈ RealDomain(; lower = 0)
        σ ∈ RealDomain(; lower = 0)
    end

    @pre begin
        CL = tvcl
        Vc = tvvc
        Ka1 = tvka1
        Ka2 = tvka2
    end

    @dosecontrol begin
        lags = (; SR = tvlag)
        bioav = (; IR = tvbio, SR = (1 - tvbio))
    end

    @dynamics begin
        IR' = -Ka1 * IR
        SR' = -Ka2 * SR
        Central' = Ka1 * IR + Ka2 * SR - Central * CL / Vc
    end

    @derived begin
        conc = @. Normal(Central / Vc, σ)
    end

end

PumasModel
  Parameters: tvcl, tvvc, tvka1, tvka2, tvlag, tvbio, σ
  Random effects: 
  Covariates: 
  Dynamical system variables: IR, SR, Central
  Dynamical system type: Matrix exponential
  Derived: conc
  Observed: conc

param4 = (; tvcl = 5, tvvc = 50, tvka1 = 0.8, tvka2 = 0.6, tvlag = 5, tvbio = 0.5, σ = 0.0)
dose_fo1 = DosageRegimen(100; cmt = 1, time = 0)
dose_fo2 = DosageRegimen(100; cmt = 2, time = 0)
dose4 = DosageRegimen(dose_fo1, dose_fo2)  # Actual dose is made up of 2 virtual doses
subj4 = map(i -> Subject(; id = i, events = dose4, observations = (; conc = nothing)), 1:10)
sims4 = simobs(two_parallel_foabs, subj4, param4; obstimes = 0:0.1:48)

Simulated population (Vector{<:Subject})
  Simulated subjects: 10
  Simulated variables: conc

sim_plot(
    sims4,
    observations = [:conc],
    figure = (; fontsize = 18),
    axis = (;
        xlabel = "Time (hr)",
        ylabel = "Predicted Concentration (ng/mL)",
        title = "Two Parallel first-order absorption",
    ),
)

6 Transit compartments

The transit compartment model simply uses the fact that if you set up a sequence of compartments with transfer rates of Ktr you get a model where the absorption of the drug is delayed. However, instead of a fixed delay like in a lag model the drug absorption is distributed across time with a peak absorption that is not at the time of dosing. For a fixed transfer rate Ktr we get an ever growing delay as the number of transit compartments increases. Hence, it is easier to interpret the models in terms of the Mean Absorption Time MAT that is defined as MAT = (n+1)/Ktr. This also means that if we use MAT as the estimated parameter, the implied Ktr for a given number of compartments is Ktr = (n+1)/MAT.

Let us look at the Pumas code for a 1-compartment model with transit compartment absorption structure.

transit3 = @model begin

    @param begin
        tvcl ∈ RealDomain(; lower = 0)
        tvvc ∈ RealDomain(; lower = 0)
        tvktr ∈ RealDomain(; lower = 0)
        σ ∈ RealDomain(; lower = 0)
    end

    @pre begin
        CL = tvcl
        Vc = tvvc
        Ktr = tvktr
    end

    @dynamics begin
        Transit1' = -Ktr * Transit1
        Transit2' = Ktr * Transit1 - Ktr * Transit2
        Transit3' = Ktr * Transit2 - Ktr * Transit3
        Central' = Ktr * Transit3 - (CL / Vc) * Central
    end

    @derived begin
        conc = @. Normal(Central / Vc, σ)
    end
end

PumasModel
  Parameters: tvcl, tvvc, tvktr, σ
  Random effects: 
  Covariates: 
  Dynamical system variables: Transit1, Transit2, Transit3, Central
  Dynamical system type: Matrix exponential
  Derived: conc
  Observed: conc

param1 = (; tvcl = 5, tvvc = 20, tvka = 1, tvktr = 3 / 2, σ = 0.0)
sims1 = simobs(transit3, subj, param1; obstimes = 0:0.1:48)

SimulatedObservations
  Simulated variables: conc
  Time: 0.0:0.1:48.0

sim_plot(
    sims1;
    observations = [:conc],
    figure = (; fontsize = 18),
    axis = (;
        xlabel = "Time (hr)",
        ylabel = "Predicted Concentration (ng/mL)",
        title = "Three transit compartment absorption",
    ),
)

The explicitly coded version of transit compartments can be instructive, but is in general not need, as we can more succinctly describe and include it using the @delay macro.

Note

Some papers like to include an Absorption compartment following the transit compartments. We do not include that here. Mainly, the Absorption compartment only adds another parameter to estimate and doesn’t add significant flexibility or realism to the model.

7 Erlang distribution

As described in [ref], the transit compartment closed form solution is equivalent to setting the absorption rate to the probability density function (pdf) of an Erlang(N, s) distribution with integer shape parameter N and real valued scale parameter s scaled by the dose amount. We can verify that by using the @delay macro. To use the @delay macro you need to

Create a dose with the amount and dose time
Pick a compartment target that is the Central compartment of interest. A Depot or Absorption compartment is generally not needed.
Pick the distribution to use for the pdf evaluations (Erlang for Transit compartments) and its parameterization.
Set the bioavailability of the dose to 0.0 for the target compartment

The parameterization of the Erlang distribution when related to the transit compartment model is that Erlang(N, s) will represent a transit model with N-1 transit compartments and a mean absorption time (MAT) of MAT = N*s. In other words, if you have a given MAT and a given number of compartments N the distribution we need is Erlang(N, MAT/N). The following panel shows an example where we formulate an equivalent model to the transit compartment model with three compartments and a Ktr of 3/2 using the @delay macro and Erlang(N, MAT/N) distribution where MAT=2 and N=3.

erlang = @model begin
    @param begin
        tvCL ∈ RealDomain(; lower = 0)
        tvVc ∈ RealDomain(; lower = 0)
        tvMAT ∈ RealDomain(; lower = 0)
        tvN ∈ RealDomain(; lower = 0)
        σ ∈ RealDomain(; lower = 0)
    end
    @pre begin
        CL = tvCL
        Vc = tvVc
        MAT = tvMAT
        N = tvN
        ρ = @delay(Erlang(N, MAT / N), Central)
    end
    @dosecontrol begin
        bioav = (; Central = 0.0)
    end
    @dynamics begin
        Central' = ρ - (CL / Vc) * Central
    end

    @derived begin
        conc = @. Normal(Central / Vc, σ)
    end
end

PumasModel
  Parameters: tvCL, tvVc, tvMAT, tvN, σ
  Random effects: 
  Covariates: 
  Dynamical system variables: Central
  Dynamical system type: Nonlinear ODE
  Derived: conc
  Observed: conc

param2 = (; tvCL = 1, tvVc = 10, tvMAT = 2, tvN = 3, σ = 0.0)
dose2 = DosageRegimen(100; time = 0)
subj2 = Subject(; id = "1", events = dose2, observations = (; conc = nothing))

sims2 = simobs(erlang, subj2, param2; obstimes = 0:0.1:48)

SimulatedObservations
  Simulated variables: conc
  Time: 0.0:0.1:48.0

sim_plot(
    erlang,
    sims2;
    observations = [:conc],
    figure = (; fontsize = 18),
    axis = (;
        xlabel = "Time (hr)",
        ylabel = "Predicted Concentration (ng/mL)",
        title = "Erlang absorption",
    ),
)

As we can see, the solution is the same as before are the same. From the example shown, it is also clear that the @delay interface is much for flexible in the number of compartments, or rather in the shape and scale of the distribution. Remember, the transit compartment structure is not inspired by physiology it is just some way of getting an absorption that is distributed across time is a less sudden manner than the lag model. This is why we emphasize N as being the shape parameter that controls how widely to spread the absorption across time and the scale parameter MAT/N that controls the mean time of any particle to reach the central compartment. Freeing ourselves from the transit compartment model opens up the possibility to look a the Gamma and Weibull models amongst others.

8 Gamma distribution

If we forget the connection to transit compartments for a moment, the Erlang distribution used above is simply the Gamma distribution where the shape parameter is restricted to be an integer. Using the Gamma distribution instead opens up the possibility to easily estimate the shape parameters as it is not a continuous variable. When N is an integer we would need to manually construct a sequence of models to look for the optimal shape or to use the more complicated Mixed Integer Programming (MIP) methods that allow you to have both continuous and integer optimization variables. Instead, if we just do not restrict ourselves to the “transit compartment” terminology, we can just use the Gamma distribution, and optimize N normally. We may get N=4.89 to be the optimal shape parameter but that doesn’t matter because the transit compartment structure was not important in itself it was only a convenient way to introduce the delay in our ODE framework.

gamma = @model begin
    @param begin
        tvCL ∈ RealDomain(; lower = 0)
        tvVc ∈ RealDomain(; lower = 0)
        tvMAT ∈ RealDomain(; lower = 0)
        tvN ∈ RealDomain(; lower = 0)
        σ ∈ RealDomain(; lower = 0)
    end
    @pre begin
        CL = tvCL
        Vc = tvVc
        MAT = tvMAT
        N = tvN
        ρ = @delay(Gamma(N, MAT / N), Central)
    end
    @dosecontrol begin
        bioav = (; Central = 0.0)
    end
    @dynamics begin
        Central' = ρ - (CL / Vc) * Central
    end

    @derived begin
        conc = @. Normal(Central / Vc, σ)
    end
end

PumasModel
  Parameters: tvCL, tvVc, tvMAT, tvN, σ
  Random effects: 
  Covariates: 
  Dynamical system variables: Central
  Dynamical system type: Nonlinear ODE
  Derived: conc
  Observed: conc

param3 = (; tvCL = 1, tvVc = 10, tvMAT = 3, tvN = 3, σ = 0.0)
dose3 = DosageRegimen(100; time = 0)
subj3 = Subject(; id = "1", events = dose3, observations = (; conc = nothing))

sims3 = simobs(gamma, subj3, param3; obstimes = 0:0.1:48)

SimulatedObservations
  Simulated variables: conc
  Time: 0.0:0.1:48.0

sim_plot(
    gamma,
    sims3;
    observations = [:conc],
    figure = (; fontsize = 18),
    axis = (;
        xlabel = "Time (hr)",
        ylabel = "Predicted Concentration (ng/mL)",
        title = "Gamma absorption",
    ),
)

You should be able to compare the plot in the panel above to the plot in the panel for the Erlang distribution and see that they are essentially the same for the N=3 specification.

9 Weibull distribution

Another popular choice is the Weibull absorption model. This is very simply implemented in Pumas by using the @delay macro and the Weibull distribution specification. This type of absorption was introduced in Pietroskii (1987). In Pumas, we parameterize the Weibull distribution with shape k and scale λ. This is in sync with the parameterization here (wikipedia). To relate these to the parameters in Piotrovskii (1987) use λₐ = (1/λ)^s and s = k. Then, if k = s = 1 we get that λₐ we actually have an exponential distribution and the absorption reduces to first order kinetics. As before, we can characterize the mean absorption time from the distribution. To output MAT we need to calculate the mean of a given specification.

weibull = @model begin
    @param begin
        tvCL ∈ RealDomain(; lower = 0)
        tvVc ∈ RealDomain(; lower = 0)
        tvk ∈ RealDomain(lower = 0.0)
        tvλ ∈ RealDomain(lower = 0.0)
        σ ∈ RealDomain(; lower = 0)
    end
    @pre begin
        CL = tvCL
        Vc = tvVc
        k = tvk
        λ = tvλ
        dist = Weibull(k, λ)
        ρ = @delay(dist, Central)
        MAT = mean(dist)
    end
    @dosecontrol begin
        bioav = (; Central = 0.0)
    end
    @dynamics begin
        Central' = ρ - (CL / Vc) * Central
    end

    @derived begin
        conc = @. Normal(Central / Vc, σ)
    end
end

PumasModel
  Parameters: tvCL, tvVc, tvk, tvλ, σ
  Random effects: 
  Covariates: 
  Dynamical system variables: Central
  Dynamical system type: Nonlinear ODE
  Derived: conc
  Observed: conc

param4 = (; tvCL = 1, tvVc = 10, tvk = 2, tvλ = 5, σ = 0.0)
dose4 = DosageRegimen(100; time = 0)
subj4 = Subject(; id = "1", events = dose4, observations = (; conc = nothing))

sims4 = simobs(weibull, subj4, param4; obstimes = 0:0.1:48)

SimulatedObservations
  Simulated variables: conc
  Time: 0.0:0.1:48.0

sim_plot(
    weibull,
    sims4;
    observations = [:conc],
    figure = (; fontsize = 18),
    axis = (;
        xlabel = "Time (hr)",
        ylabel = "Predicted Concentration (ng/mL)",
        title = "Weibull absorption",
    ),
)

first(DataFrame(sims4), 20)[!, [:time, :conc, :MAT]]

20×3 DataFrame

Row	time	conc	MAT
	Float64	Float64?	Float64?
1	0.0	missing	4.43113
2	0.0	0.0	4.43113
3	0.1	0.0039859	4.43113
4	0.2	0.0158811	4.43113
5	0.3	0.0355783	4.43113
6	0.4	0.0629525	4.43113
7	0.5	0.0978607	4.43113
8	0.6	0.140143	4.43113
9	0.7	0.189624	4.43113
10	0.8	0.24611	4.43113
11	0.9	0.309396	4.43113
12	1.0	0.379259	4.43113
13	1.1	0.455461	4.43113
14	1.2	0.53775	4.43113
15	1.3	0.625865	4.43113
16	1.4	0.719536	4.43113
17	1.5	0.818481	4.43113
18	1.6	0.922408	4.43113
19	1.7	1.03102	4.43113
20	1.8	1.144	4.43113

10 Conclusion

We have seen in this tutorial some examples of how to implement various delayed absorption models that are commonly employed in pharmacokinetic modeling using Pumas. The powerful Pumas modeling language framework makes it relatively straightforward to implement these and other complex PK and PKPD models for modeling and simulation purposes.

This concludes our tutorial on advanced absorption models with Pumas. Thanks for reading, and please be sure to check out our other tutorials as well.

Reuse

CC BY-SA 4.0