using Random
using Pumas
using PumasUtilities
using CairoMakie
using AlgebraOfGraphics
using DataFramesMeta
using CSV
using Dates
PK37 - Invitro/Invivo Extrapolation I
1 Learning Outcome
In this exercise, you will learn to assess the relationship between rate and concentration using a non-linear clearance model. The Clearance is impacted by two parallel non-linear processes, and both of them are taken into account for calculating rates. The use of a single enzyme system was not able to account for the deviations in the system.
2 Libraries
Call the necessary libraries to get started
3 Model
pk_37 = @model begin
@metadata begin
desc = "Invivo/Invitro Extrapolation Model"
end
@param begin
"""
Maximum Rate of System 1 (μmol/min/g)
"""
tvvmax1 ∈ RealDomain(lower = 0)
"""
Michaelis-Menten constant for system 1 (μmol/L)
"""
tvkm1 ∈ RealDomain(lower = 0)
"""
Maximum Rate of System 2 (μmol/min/g)
"""
tvvmax2 ∈ RealDomain(lower = 0)
"""
Michaelis-Menten constant for system 2 (μmol/L)
"""
tvkm2 ∈ RealDomain(lower = 0)
Ω ∈ PDiagDomain(4)
"""
Additive RUV
"""
σ_add ∈ RealDomain(lower = 0)
end
@random begin
η ~ MvNormal(Ω)
end
@covariates Conc
@pre begin
Vmax1 = tvvmax1 * exp(η[1])
Km1 = tvkm1 * exp(η[2])
Vmax2 = tvvmax2 * exp(η[3])
Km2 = tvkm2 * exp(η[4])
_Conc = Conc
end
@derived begin
rate = @. ((Vmax1 * _Conc / (Km1 + _Conc)) + (Vmax2 * _Conc / (Km2 + _Conc)))
"""
Observed Concentrations (μmol/L)
"""
dv_rate ~ @. Normal(rate, σ_add)
end
endPumasModel
Parameters: tvvmax1, tvkm1, tvvmax2, tvkm2, Ω, σ_add
Random effects: η
Covariates: Conc
Dynamical system variables:
Dynamical system type: No dynamical model
Derived: rate, dv_rate
Observed: rate, dv_rate
4 Parameters
The parameters are as given below. Note that tv represents the typical value for parameters.
tvvmax1- Maximum Rate of System 1 (μmol/min/g)tvkm1- Michaelis-Menten constant for system 1 (μmol/L)tvvmax2- Maximum Rate of System 2 (μmol/min/g)tvkm2- Michaelis-Menten constant for system 2 (μmol/L)
param = (
tvvmax1 = 2.43386,
tvkm1 = 256.384,
tvvmax2 = 0.22523,
tvkm2 = 2.23625,
Ω = Diagonal([0.04, 0.04, 0.04, 0.04]),
σ_add = 0.0441543,
)5 Dosage Regimen
Since there is no dosage regimen, we will create a DataFrame with no events and read the file using read_pumas which will be used for simulation.
df_sub1 = DataFrame(id = 1, time = 0:1:40000, Conc = 0.00:0.05:2000, dv_rate = missing);sub1 =
read_pumas(df_sub1, observations = [:dv_rate], covariates = [:Conc], event_data = false)Population
Subjects: 1
Covariates: Conc
Observations: dv_rate
6 Simulations
We will now simulate the rate using the simobs function
Random.seed!(123)The random effects are zero’ed out since we are simulating a single subject
zfx = zero_randeffs(pk_37, sub1, param)1-element Vector{NamedTuple{(:η,), Tuple{Vector{Float64}}}}:
(η = [0.0, 0.0, 0.0, 0.0],)
sim_sub1 = simobs(pk_37, sub1, param, zfx)Simulated population (Vector{<:Subject})
Simulated subjects: 1
Simulated variables: rate, dv_rate
df37 = DataFrame(sim_sub1)
first(df37, 5)| Row | id | time | rate | dv_rate | evid | Conc | η_1 | η_2 | η_3 | η_4 | Vmax1 | Km1 | Vmax2 | Km2 | _Conc |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| String? | Float64 | Float64? | Float64? | Int64? | Float64? | Float64? | Float64? | Float64? | Float64? | Float64? | Float64? | Float64? | Float64? | Float64? | |
| 1 | 1 | 0.0 | 0.0 | -0.000211197 | 0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 2.43386 | 256.384 | 0.22523 | 2.23625 | 0.0 |
| 2 | 1 | 1.0 | 0.00540031 | -0.00792422 | 0 | 0.05 | 0.0 | 0.0 | 0.0 | 0.0 | 2.43386 | 256.384 | 0.22523 | 2.23625 | 0.05 |
| 3 | 1 | 2.0 | 0.0105896 | -0.0176694 | 0 | 0.1 | 0.0 | 0.0 | 0.0 | 0.0 | 2.43386 | 256.384 | 0.22523 | 2.23625 | 0.1 |
| 4 | 1 | 3.0 | 0.0155811 | 0.00176718 | 0 | 0.15 | 0.0 | 0.0 | 0.0 | 0.0 | 2.43386 | 256.384 | 0.22523 | 2.23625 | 0.15 |
| 5 | 1 | 4.0 | 0.020387 | 0.0739383 | 0 | 0.2 | 0.0 | 0.0 | 0.0 | 0.0 | 2.43386 | 256.384 | 0.22523 | 2.23625 | 0.2 |
7 Visualization
Use the DataFrame for plotting
@chain df37 begin
@rsubset :Conc ∈
[0.02, 0.05, 0.1, 0.2, 0.5, 1, 10, 20, 50, 100, 180, 200, 250, 500, 1000, 2000]
data(_) *
mapping(:Conc => "Concentration (μM)", :rate => "Rate (μmol/g/min)") *
visual(ScatterLines; linewidth = 4, markersize = 12, markercolor = :blue)
draw(
axis = (;
title = "Metabolic Rate vs Substrate Concentration",
xticks = 0:200:2000,
yticks = 0:0.5:2.5,
),
)
end@chain df37 begin
@rsubset :Conc ∈
[0.02, 0.05, 0.1, 0.2, 0.5, 1, 10, 20, 50, 100, 180, 200, 250, 500, 1000, 2000]
data(_) *
mapping(:Conc => "Concentration (μM)", :rate => "Rate (μmol/g/min)") *
visual(ScatterLines; linewidth = 4, markersize = 12, markercolor = :blue)
draw(
axis = (;
title = "Metabolic Rate vs Substrate Concentration",
yticks = 0:0.5:2.5,
xscale = log10,
),
)
end8 Population simulation
par = (
tvvmax1 = 2.43386,
tvkm1 = 256.384,
tvvmax2 = 0.22523,
tvkm2 = 2.23625,
Ω = Diagonal([0.0426, 0.0124, 0.0256, 0.0371]),
σ_add = 0.0441543,
)
df1_pop = map(
i -> DataFrame(
id = i,
time = (1:1:16),
Conc = [
0.02,
0.05,
0.1,
0.2,
0.5,
1,
10,
20,
50,
100,
180,
200,
250,
500,
1000,
2000,
],
dv_rate = missing,
),
1:55,
)
df2_pop = vcat(DataFrame.(df1_pop)...)
pop =
read_pumas(df2_pop, observations = [:dv_rate], covariates = [:Conc], event_data = false)
Random.seed!(1234)
sim_pop = simobs(pk_37, pop, par)
df_sim = DataFrame(sim_pop);
#CSV.write("pk_37.csv", df_sim)