using Random
using Pumas
using PumasUtilities
using CairoMakie
using AlgebraOfGraphics
using CSV
using DataFramesMeta
using Dates
PK31 - Turnover II - Intravenous Dosing of Hormone
1 Background
- Structural model - Two compartment with additional input for basal hormone synthesis in the central compartment
- Route of administration - IV infusion (1 minute)
- Dosage Regimen - 36,630 pmol
- Number of Subjects - 1
2 Learning Outcome
In this model, you will learn how to build a two compartment model with additional input for basal hormone level. This model will help to simulate the plasma concentration profile after IV administration, considering basal hormone input.
3 Objectives
- To analyze the intravenous datasets with parallel turnover
- To write a multi-compartment model in terms of differential equations
4 Libraries
Call the necessary libraries to get start.
5 Model
In this one compartment model, we administer an IV infusion dose to the central compartment.
= @model begin
pk_31 @metadata begin
= "Two Compartment Model"
desc = u"hr"
timeu end
@param begin
"""
Basal hormonal input (pmol/hr)
"""
∈ RealDomain(lower = 0)
tvkin """
Central Volume of Distribution (L)
"""
∈ RealDomain(lower = 0)
tvvc """
Clearance (L/hr)
"""
∈ RealDomain(lower = 0)
tvcl """
Intercompartmental Clearance (L/hr)
"""
∈ RealDomain(lower = 0)
tvq """
Peripheral Volume of Distribution (L)
"""
∈ RealDomain(lower = 0)
tvvp ∈ PDiagDomain(5)
Ω """
Proportional RUV
"""
∈ RealDomain(lower = 0)
σ²_prop
end
@random begin
~ MvNormal(Ω)
η end
@pre begin
= tvkin * exp(η[1])
Kin = tvvc * exp(η[2])
Vc = tvcl * exp(η[3])
Cl = tvq * exp(η[4])
Q = tvvp * exp(η[5])
Vp end
@dynamics begin
' = Kin - (Cl / Vc) * Central - (Q / Vc) * Central + (Q / Vp) * Peripheral
Central' = (Q / Vc) * Central - (Q / Vp) * Peripheral
Peripheralend
@derived begin
= @. Central / Vc
cp """
Observed Concentration (pmol/L)
"""
~ @. Normal(cp, sqrt(cp^2 * σ²_prop))
dv end
end
PumasModel
Parameters: tvkin, tvvc, tvcl, tvq, tvvp, Ω, σ²_prop
Random effects: η
Covariates:
Dynamical system variables: Central, Peripheral
Dynamical system type: Matrix exponential
Derived: cp, dv
Observed: cp, dv
6 Parameters
The parameters are as given below. Note that tv
represents the typical value for parameters.
Kin
- Basal hormonal input (pmol/hr)Vc
- Central Volume of Distribution (L)Cl
- Clearance (L/hr)Q
- Intercompartmental Clearance (L/hr)Vp
- Peripheral Volume of Distribution (L)Ω
- Between Subject Variability,σ
- Residual error
= (
param = 1531.87,
tvkin = 8.8455,
tvvc = 76.5987,
tvcl = 56.8775,
tvq = 58.8033,
tvvp = Diagonal([0.04, 0.04, 0.04, 0.04, 0.04]),
Ω = 0.015,
σ²_prop )
7 Dosage Regimen
A single dose of 36630 pmol is given as an IV Infusion
= DosageRegimen(36630, time = 0, cmt = 1, duration = 0.0166) ev1
1×10 DataFrame
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 | 36630.0 | 1 | 0.0 | 0 | 2.20663e6 | 0.0166 | 0 | NullRoute |
= Subject(id = 1, events = ev1, observations = (cp = nothing,)) sub1
Subject
ID: 1
Events: 1
Observations: cp: (n=0)
8 Simulation
Simulate the plasma concentration profile
Random.seed!(123)
The random effects are zero’ed out since we are simulating a single subject
= zero_randeffs(pk_31, sub1, param) zfx
(η = [0.0, 0.0, 0.0, 0.0, 0.0],)
= simobs(pk_31, sub1, param, zfx, obstimes = 0.01:0.0001:32) sim_sub1
SimulatedObservations
Simulated variables: cp, dv
Time: 0.01:0.0001:32.0
9 Visualization
@chain DataFrame(sim_sub1) begin
dropmissing(:cp)
data(_) *
mapping(:time => "Time (hours)", :cp => "Concentration (pmol/L)") *
visual(Lines; linewidth = 4)
draw(;
= (; fontsize = 22),
figure = (;
axis = log10,
yscale = i -> (@. string(round(i; digits = 1))),
ytickformat = 0:5:35,
xticks
),
)end
10 Population simulation
= (
par = 1531.87,
tvkin = 8.8455,
tvvc = 76.5987,
tvcl = 56.8775,
tvq = 58.8033,
tvvp = Diagonal([0.09, 0.0125, 0.0225, 0.04, 0.0365]),
Ω = 0.0612144,
σ²_prop
)
= DosageRegimen(36630, time = 0, cmt = 1, duration = 0.0166)
ev1 = map(i -> Subject(id = i, events = ev1), 1:85)
pop
= simobs(
sim_pop
pk_31,
pop,
par,= [
obstimes 0.0167,
0.1167,
0.167,
0.25,
0.583,
0.833,
1.083,
1.583,
2.083,
4.083,
8.083,
12,
23.5,
24.25,
26.75,
32,
],
)
= DataFrame(sim_pop)
df_sim #CSV.write("pk_31.csv", df_sim)