```
using Random
using Pumas
using PumasUtilities
using CairoMakie
Random.seed!(1234) # Set random seed for reproducibility
```

# Absorption Models

## 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 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:

```
= @model begin
foabs
@param begin
∈ RealDomain(; lower = 0)
tvcl ∈ RealDomain(; lower = 0)
tvvc ∈ RealDomain(; lower = 0)
tvka ∈ RealDomain(; lower = 0)
σ end
@pre begin
= tvcl
CL = tvvc
Vc = tvka
Ka end
@dynamics begin
' = -Ka * Depot
Depot' = Ka * Depot - (CL / Vc) * Central
Centralend
@derived begin
= @. Normal(Central / Vc, σ)
conc end
end
```

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

```
= (; tvcl = 5, tvvc = 20, tvka = 1, σ = 0.0)
param1
= DosageRegimen(100; time = 0)
dose1 = Subject(; id = "1", events = dose1, observations = (; conc = nothing))
subj
= simobs(foabs, subj, param1; obstimes = 0:0.1:48) sims1
```

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

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

## 3 Zero-Order Absorption

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:

```
= @model begin
zoabs @param begin
∈ RealDomain(; lower = 0)
tvcl ∈ RealDomain(; lower = 0)
tvvc ∈ RealDomain(; lower = 0)
tvdur ∈ RealDomain(; lower = 0)
σ end
@pre begin
= tvcl
CL = tvvc
Vc end
@dosecontrol begin
= (; Central = tvdur)
duration end
@dynamics begin
' = -(CL / Vc) * Central
Centralend
@derived begin
= @. Normal(Central / Vc, σ)
conc end
end
```

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

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

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

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

## 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).

```
= @model begin
zofo_paral_abs
@param begin
∈ RealDomain(; lower = 0)
tvcl ∈ RealDomain(; lower = 0)
tvvc ∈ RealDomain(; lower = 0)
tvka ∈ RealDomain(; lower = 0)
tvdur ∈ RealDomain(; lower = 0)
tvbio ∈ RealDomain(; lower = 0)
tvlag ∈ RealDomain(; lower = 0)
σ end
@pre begin
= tvcl
CL = tvvc
Vc = tvka
Ka end
@dosecontrol begin
= (; Central = tvdur)
duration = (; Depot = tvbio, Central = 1 - tvbio)
bioav = (; Depot = tvlag)
lags end
@dynamics begin
' = -Ka * Depot
Depot' = Ka * Depot - (CL / Vc) * Central
Centralend
@derived begin
= @. Normal(Central / Vc, σ)
conc 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
```

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.

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

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

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

## 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:

```
= @model begin
two_parallel_foabs
@param begin
∈ RealDomain(; lower = 0)
tvcl ∈ RealDomain(; lower = 0)
tvvc ∈ RealDomain(; lower = 0)
tvka1 ∈ RealDomain(; lower = 0)
tvka2 ∈ RealDomain(; lower = 0)
tvlag ∈ RealDomain(; lower = 0)
tvbio ∈ RealDomain(; lower = 0)
σ end
@pre begin
= tvcl
CL = tvvc
Vc = tvka1
Ka1 = tvka2
Ka2 end
@dosecontrol begin
= (; SR = tvlag)
lags = (; IR = tvbio, SR = (1 - tvbio))
bioav end
@dynamics begin
' = -Ka1 * IR
IR' = -Ka2 * SR
SR' = Ka1 * IR + Ka2 * SR - Central * CL / Vc
Centralend
@derived begin
= @. Normal(Central / Vc, σ)
conc 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
```

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

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

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

## 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.

```
= @model begin
transit3
@param begin
∈ RealDomain(; lower = 0)
tvcl ∈ RealDomain(; lower = 0)
tvvc ∈ RealDomain(; lower = 0)
tvktr ∈ RealDomain(; lower = 0)
σ end
@pre begin
= tvcl
CL = tvvc
Vc = tvktr
Ktr end
@dynamics begin
' = -Ktr * Transit1
Transit1' = Ktr * Transit1 - Ktr * Transit2
Transit2' = Ktr * Transit2 - Ktr * Transit3
Transit3' = Ktr * Transit3 - (CL / Vc) * Central
Centralend
@derived begin
= @. Normal(Central / Vc, σ)
conc 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
```

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

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

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

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.

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`

.

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

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

```
= (; tvCL = 1, tvVc = 10, tvMAT = 2, tvN = 3, σ = 0.0)
param2 = DosageRegimen(100; time = 0)
dose2 = Subject(; id = "1", events = dose2, observations = (; conc = nothing))
subj2
= simobs(erlang, subj2, param2; obstimes = 0:0.1:48) sims2
```

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

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

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.

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

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

```
= (; tvCL = 1, tvVc = 10, tvMAT = 3, tvN = 3, σ = 0.0)
param3 = DosageRegimen(100; time = 0)
dose3 = Subject(; id = "1", events = dose3, observations = (; conc = nothing))
subj3
= simobs(gamma, subj3, param3; obstimes = 0:0.1:48) sims3
```

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

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

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.

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

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

```
= (; tvCL = 1, tvVc = 10, tvk = 2, tvλ = 5, σ = 0.0)
param4 = DosageRegimen(100; time = 0)
dose4 = Subject(; id = "1", events = dose4, observations = (; conc = nothing))
subj4
= simobs(weibull, subj4, param4; obstimes = 0:0.1:48) sims4
```

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

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

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

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.