using Random
using Pumas
using PumasUtilities
using AlgebraOfGraphics
using CairoMakie
using CSV
using DataFramesMeta
using Dates
PK16 - Two-compartment intravenous plasma/urine
1 Background
- Structural Model - Two Compartment model with first order elimination
- Route of Administration - Intravenous infusion (Multiple-dose)
- Dosage Regimen - 538 µmol/kg followed by 3390 µmol/kg
- Number of Subjects - 1
2 Learning Outcome
This exercise explains a simultaneous analysis of plasma and urine data using a two-compartment model, with an additional urinary compartment accounting for the amount of drug excreted in urine.
3 Objectives
In this exercise, you will learn how to build a two compartment model and to simulate the data for a single subject based on the given dosage regimen.
4 Libraries
Call the necessary libraries to get started
5 Model
A two-compartment model with one Central and one Peripheral compartment will help in understanding the plasma concentration. A separate Urine compartment accounts for the fraction (amount) of drug excreted in urine.
= @model begin
pk_16 @metadata begin
= "Two Compartment Model with Urine Compartment"
desc = u"hr"
timeu end
@param begin
"""
Renal Clearance (L/hr/kg)
"""
∈ RealDomain(lower = 0)
tvclr """
Non-renal Clearance (L/hr/kg)
"""
∈ RealDomain(lower = 0)
tvclm """
Volume of Central Compartment (L/kg)
"""
∈ RealDomain(lower = 0)
tvvc """
Volume of Peripheral Compartment (L/kg)
"""
∈ RealDomain(lower = 0)
tvvp """
Inter-compartmental Clearance (L/hr/kg)
"""
∈ RealDomain(lower = 0)
tvq ∈ PDiagDomain(5)
Ω """
Proportional RUV - Plasma
"""
∈ RealDomain(lower = 0)
σ₁_prop """
Proportional RUV - Urine
"""
∈ RealDomain(lower = 0)
σ₂_prop end
@random begin
~ MvNormal(Ω)
η end
@pre begin
= tvclr * exp(η[1])
CLr = tvclm * exp(η[2])
CLm = tvvc * exp(η[3])
Vc = tvvp * exp(η[4])
Vp = tvq * exp(η[5])
Q = CLr + CLm
CL = Vc + Vp
Vss = 0.693 * Vss / CL
t12 end
@dynamics begin
' =
Central-(CLr / Vc) * Central - (CLm / Vc) * Central + (Q / Vp) * Peripheral -
/ Vc) * Central
(Q ' = (Q / Vc) * Central - (Q / Vp) * Peripheral
Peripheral' = (CLr / Vc) * Central
Urineend
@derived begin
= @. Central / (Vc + eps())
cp_plasma """
Observed Plasma Concentration (μM)
"""
~ @. Normal(cp_plasma, abs(cp_plasma) * σ₁_prop)
dv_plasma = @. Urine
cp_urine """
Observed Urine Amount (μmol)
"""
~ @. Normal(cp_urine, abs(cp_urine) * σ₂_prop)
dv_urine end
end
PumasModel
Parameters: tvclr, tvclm, tvvc, tvvp, tvq, Ω, σ₁_prop, σ₂_prop
Random effects: η
Covariates:
Dynamical system variables: Central, Peripheral, Urine
Dynamical system type: Matrix exponential
Derived: cp_plasma, dv_plasma, cp_urine, dv_urine
Observed: cp_plasma, dv_plasma, cp_urine, dv_urine
6 Parameters
CLr
- Renal Clearance (L/hr/kg)CLm
- Non-renal Clearance (L/hr/kg)Vc
- Volume of Central Compartment (L/kg)Vp
- Volume of Peripheral Compartment (L/kg)Q
- Inter-compartmental Clearance (L/hr/kg)Ω
- Between Subject Variabilityσ²₁
- Residual error 1 (for plasma conc)σ²₂
- Residual error 2 (for amount excreted in urine)
Derived Parameters:
CL
- Total Clearance (L/hr/kg)Vss
- Volume of distribution at steady state (L/kg)t1/2
- Half life (hr)
Typical value (tv
) estimates for individual parameters
= (
param = 0.05,
tvclm = 0.31,
tvclr = 1.6,
tvvc = 0.16,
tvvp = 0.03,
tvq = Diagonal([0.04, 0.04, 0.04, 0.04, 0.04]),
Ω = 0.1,
σ₁_prop = 0.1,
σ₂_prop )
7 Dosage Regimen
- The subject (male dog) received two consecutive intravenous infusions of a drug.
- The subject received an initial dose of 538 µmol/kg from time 0 to 0.983 hr (rate = 547.304) followed by 3390 µmol/kg dose from time 0.983 to 23.95 hr (rate = 147.603).
- Total infused dose: 3928 µmol/kg.
= DosageRegimen([538, 3390], time = [0, 0.983], cmt = [1, 1], rate = [547.304, 147.603]) ev1
Row | time | cmt | amt | evid | ii | addl | rate | duration | ss | route |
---|---|---|---|---|---|---|---|---|---|---|
Float64 | Int64 | Float64 | Int8 | Float64 | Int64 | Float64 | Float64 | Int8 | NCA.Route | |
1 | 0.0 | 1 | 538.0 | 1 | 0.0 | 0 | 547.304 | 0.983 | 0 | NullRoute |
2 | 0.983 | 1 | 3390.0 | 1 | 0.0 | 0 | 147.603 | 22.967 | 0 | NullRoute |
= Subject(id = 1, events = ev1, observations = (cp = nothing,)) sub1
Subject
ID: 1
Events: 4
Observations: cp: (n=0)
8 Simulation
To simulate the following for the above subject with specific observation time points.
- Plasma concentration.
- Amount excreted unchanged in urine.
Random.seed!(123)
The random effects are zero’ed out since we are simulating means
= zero_randeffs(pk_16, sub1, param) zfx
(η = [0.0, 0.0, 0.0, 0.0, 0.0],)
= simobs(
sim_sub1
pk_16,
sub1,
param,
zfx,= [
obstimes 0.5,
1,
2,
4,
6.1,
7.6,
8.02,
12.05,
12.15,
15.95,
22.13,
23.89,
24.05,
24.46,
24.94,
25.94,
26.96,
27.95,
29.97,
31.94,
35.96,
36.2,
48,
48.2,
54,
60,
60.2,
72,
72.2,
], )
SimulatedObservations
Simulated variables: cp_plasma, dv_plasma, cp_urine, dv_urine
Time: [0.5, 1.0, 2.0, 4.0, 6.1, 7.6, 8.02, 12.05, 12.15, 15.95 … 31.94, 35.96, 36.2, 48.0, 48.2, 54.0, 60.0, 60.2, 72.0, 72.2]
9 Visualization
Create a DataFrame
of the simulated data.
= DataFrame(sim_sub1)
df1 first(df1, 5)
Row | id | time | cp_plasma | dv_plasma | cp_urine | dv_urine | evid | amt | cmt | rate | duration | ss | ii | route | η_1 | η_2 | η_3 | η_4 | η_5 | Central | Peripheral | Urine | CLr | CLm | Vc | Vp | Q | CL | Vss | t12 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
String? | Float64 | Float64? | Float64? | Float64? | Float64? | Int64? | Float64? | Symbol? | Float64? | Float64? | Int8? | Float64? | NCA.Route? | Float64? | Float64? | Float64? | Float64? | Float64? | Float64? | Float64? | Float64? | Float64? | Float64? | Float64? | Float64? | Float64? | Float64? | Float64? | Float64? | |
1 | 1 | 0.0 | missing | missing | missing | missing | 1 | 538.0 | Central | 547.304 | 0.983 | 0 | 0.0 | NullRoute | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | missing | missing | missing | 0.31 | 0.05 | 1.6 | 0.16 | 0.03 | 0.36 | 1.76 | 3.388 |
2 | 1 | 0.5 | 161.044 | 166.797 | 12.7335 | 11.7625 | 0 | missing | missing | missing | missing | missing | missing | missing | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 257.67 | 1.19427 | 12.7335 | 0.31 | 0.05 | 1.6 | 0.16 | 0.03 | 0.36 | 1.76 | 3.388 |
3 | 1 | 0.983 | missing | missing | missing | missing | 1 | 3390.0 | Central | 147.603 | 22.967 | 0 | 0.0 | NullRoute | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | missing | missing | missing | 0.31 | 0.05 | 1.6 | 0.16 | 0.03 | 0.36 | 1.76 | 3.388 |
4 | 1 | 1.0 | 299.498 | 322.168 | 48.9648 | 39.3479 | 0 | missing | missing | missing | missing | missing | missing | missing | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 479.197 | 4.45023 | 48.9648 | 0.31 | 0.05 | 1.6 | 0.16 | 0.03 | 0.36 | 1.76 | 3.388 |
5 | 1 | 2.0 | 317.468 | 271.27 | 144.684 | 135.949 | 0 | missing | missing | missing | missing | missing | missing | missing | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 507.948 | 12.1436 | 144.684 | 0.31 | 0.05 | 1.6 | 0.16 | 0.03 | 0.36 | 1.76 | 3.388 |
= [
time_plasma 0.5,
1,
2,
4,
7.6,
8.02,
12.05,
15.95,
22.13,
23.89,
24.46,
24.94,
25.94,
26.96,
27.95,
29.97,
31.94,
35.96,
48,
54,
60,
72,
]
= @rsubset df1 :time in time_plasma df_plasma
= [6.1, 12.15, 24.05, 36.2, 48.2, 60.2, 72.2]
time_urine = @rsubset df1 :time in time_urine df_urine
Plot the simulated plasma concentration data and the amount of drug excreted in urine in a single plot.
= @chain df_plasma begin
plasma_plt @rtransform :Matrix = "Plasma"
data(_) *
mapping(
:time => "Time (hours)",
:cp_plasma => "Concentration (μM) & Amount (μmol)",
= :Matrix,
color *
) visual(ScatterLines, color = :blue, linewidth = 4)
end
= @chain df_urine begin
urine_plt @rtransform :Matrix = "Urine"
data(_) *
mapping(
:time => "Time (hours)",
:cp_urine => "Concentration (μM) & Amount (μmol)",
= :Matrix,
color *
) visual(ScatterLines, color = :black, linewidth = 4)
end
draw(
+ urine_plt;
plasma_plt = (;
axis = [0, 10, 20, 30, 40, 50, 60, 70, 80],
xticks = [0.1, 0, 1, 10, 100, 1000, 10000],
yticks = log,
yscale = x -> string.(round.(x; digits = 1)),
ytickformat = 3,
ygridwidth = true,
yminorticksvisible = true,
yminorgridvisible = IntervalsBetween(10),
yminorticks = true,
xminorticksvisible = true,
xminorgridvisible = IntervalsBetween(5),
xminorticks
),= (; fontsize = 22,),
figure )
10 Population simulation
= (
par = 0.05,
tvclm = 0.31,
tvclr = 1.6,
tvvc = 0.16,
tvvp = 0.03,
tvq = Diagonal([0.0081, 0.004, 0.0004, 0.0256, 0.0676]),
Ω = 0.04,
σ₁_prop = 0.09,
σ₂_prop
)
= DosageRegimen([538, 3390], time = [0, 0.983], cmt = [1, 1], rate = [547.304, 147.603])
ev1 = map(i -> Subject(id = i, events = ev1), 1:80)
pop
Random.seed!(1234)
= simobs(
pop_sim
pk_16,
pop,
par,= [
obstimes 0.5,
1,
2,
4,
6.1,
7.6,
8.02,
12.05,
15.95,
22.13,
23.89,
24.46,
24.94,
25.94,
26.96,
27.95,
29.97,
31.94,
36,
48,
54,
60,
72,
],
)
= DataFrame(pop_sim)
df_sim #CSV.write("pk_16.csv", df_sim)