using Random
using Pumas
using PumasUtilities
using AlgebraOfGraphics
using CairoMakie
using CSV
using DataFramesMeta
using Dates
PK24 - Non-linear Kinetics - Flow II
1 Background
- Structural model - Multi-compartment (three) model with concentration dependent clearance
- Route of administration - IV infusion
- Dosage Regimen - One IV infusion dose of 10 mg/kg given for 2 hours.
- Number of Subjects - 1
2 Learning Outcomes
This exercise explains a multi-compartment model following concentration dependent clearance.
3 Objectives
In this tutorial, you will learn:
- To build a multi-compartment model with concentration dependent clearance.
- To simulate the data for one subject given an IV infusion for 2 hours.
4 Libraries
Call the necessary libraries to get started.
5 Model
This model contains three compartments (Central, shallow and deep) and the clearance is dependent on the change in plasma concentration due to the drug, which is known to reduce cardiac output and hepatic blood flow with an increase in plasma concentration. The change in clearance can be accounted by the equation: \(CL= CLo-A(Central/Vc)\) where \(A\) is proportionality constant between \(CL\) and \(C\)
= @model begin
pk_24 @metadata begin
= "Three Compartment Model"
desc = u"hr"
timeu end
@param begin
"""
Volume of Central Compartment(L/kg)
"""
∈ RealDomain(lower = 0)
tvvc """
Clearance(L/hr/kg)
"""
∈ RealDomain(lower = 0)
tvclo """
Inter-compartmental (Shallow) Clearance(L/hr/kg)
"""
∈ RealDomain(lower = 0)
tvq1 """
Inter-compartmental (Deep) Clearance(L/hr/kg)
"""
∈ RealDomain(lower = 0)
tvq2 """
Volume of Shallow Compartment(L/kg)
"""
∈ RealDomain(lower = 0)
tvvp1 """
Volume of Deep Compartment(L/kg)
"""
∈ RealDomain(lower = 0)
tvvp2 """
Proportionality constant Between Cp and CL(L2/hr/μg/kg)
"""
∈ RealDomain(lower = 0)
tva ∈ PDiagDomain(7)
Ω ∈ RealDomain(lower = 0)
σ_prop end
@random begin
~ MvNormal(Ω)
η end
@pre begin
= tvvc * exp(η[1])
Vc = tvclo * exp(η[2])
CLo = tvq1 * exp(η[3])
Q1 = tvq2 * exp(η[4])
Q2 = tvvp1 * exp(η[5])
Vp1 = tvvp2 * exp(η[6])
Vp2 = 1000 * tva * exp(η[7])
A end
@dynamics begin
' =
Central-(CLo - (A * (Central / Vc))) * (Central / Vc) - Q1 * (Central / Vc) +
* (Shallow / Vp1) - Q2 * (Central / Vc) + Q2 * (Deep / Vp2)
Q1 ' = Q1 * (Central / Vc) - Q1 * (Shallow / Vp1)
Shallow' = Q2 * (Central / Vc) - Q2 * (Deep / Vp2)
Deepend
@derived begin
= @. 1000 * Central / Vc
cp """
Observed Concentration (μg/L)
"""
~ @. Normal(cp, sqrt(cp^2 * σ_prop))
dv end
end
PumasModel
Parameters: tvvc, tvclo, tvq1, tvq2, tvvp1, tvvp2, tva, Ω, σ_prop
Random effects: η
Covariates:
Dynamical system variables: Central, Shallow, Deep
Dynamical system type: Nonlinear ODE
Derived: cp, dv
Observed: cp, dv
6 Parameters
The parameters are as given below. Note that tv
represents the typical value for parameters.
Vc
- Volume of Central Compartment (L/kg)CLo
- Clearance (L/hr/kg)Q1
- Inter-compartmental (Shallow) Clearance (L/hr/kg)Q2
- Inter-compartmental (Deep) Clearance (L/hr/kg)Vp1
- Volume of Shallow Compartment (L/kg)Vp2
- Volume of Deep Compartment (L/kg)A
- proportionality constant Between Cp and CL (L2/hr/μg/kg)Ω
- Between Subject Variabilityσ
- Residual error (for plasma concentration)
= (
param = 0.68,
tvvc = 6.61,
tvclo = 5.94,
tvq1 = 0.93,
tvq2 = 1.77,
tvvp1 = 3.20,
tvvp2 = 0.0025,
tva = Diagonal([0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04]),
Ω = 0.004,
σ_prop )
7 Dosage Regimen
Build the dosage regimen for one subject given an IV infusion of 10 mg/kg for 2 hours.
= DosageRegimen(10, time = 0, cmt = 1, duration = 2) DR
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 | 10.0 | 1 | 0.0 | 0 | 5.0 | 2.0 | 0 | NullRoute |
= Subject(id = 1, events = DR) s1
Subject
ID: 1
Events: 2
8 Simulation
The random effects are zero’ed out since we are simulating means
= zero_randeffs(pk_24, s1, param) zfx
(η = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],)
To simulate the data with specific data points.
= simobs(pk_24, s1, param, zfx, obstimes = 0.1:0.01:8) sim_sub1
SimulatedObservations
Simulated variables: cp, dv
Time: 0.1:0.01:8.0
9 Visualization
@chain DataFrame(sim_sub1) begin
dropmissing(:cp)
data(_) *
mapping(:time => "Time (hour)", :cp => "PK24 Concentrations (μg/mL)") *
visual(Lines, linewidth = 4)
draw(; axis = (; xticks = 0:2:8), figure = (; fontsize = 22))
end
10 Population simulation
= (
parameters = 0.68,
tvvc = 6.61,
tvclo = 5.94,
tvq1 = 0.93,
tvq2 = 1.77,
tvvp1 = 3.20,
tvvp2 = 0.0025,
tva = Diagonal([0.02, 0.0, 0.02, 0.04, 0.02, 0.02, 0.0]),
Ω = 0.039,
σ_prop
)
= DosageRegimen(10, time = 0, cmt = 1, duration = 2)
DR = map(i -> Subject(id = i, events = DR), 1:45)
pop
Random.seed!(1234)
= simobs(
ss
pk_24,
pop,
parameters,= [
obstimes 0.25,
0.5,
0.75,
1.05,
1.25,
1.49,
1.75,
1.99,
2.16,
2.35,
2.4,
2.65,
2.81,
2.95,
3.11,
3.56,
4.15,
6,
7,
],
)
= DataFrame(ss)
df_sim #CSV.write("pk_24.csv", df_sim)