Commit 0e4d234a authored by Amina Gojak's avatar Amina Gojak
Browse files

Delete eleveld.jl

parent 55272818
using Turing, Distributions, DifferentialEquations, MCMCChains, Plots, StatsPlots
# Set a seed for reproducibility.
using Random
Random.seed!(14);
using Logging
Logging.disable_logging(Logging.Warn)
LogLevel(1001)
# individual values (female, 40 yo, 168 cm, 60 kg)
AGE = 40
PMA = 40+40*52 # not born prematurely and now 40 yo
weeks_40 = 0
WGT = 60
BMI = 24.1
# ref values from reference person (male, 35 yo, 170 cm, 70 kg)
AGE_ref = 35
WGT_ref = 70
PMA_ref = 40+35*52 # not born prematurely and now 35 yo
BMI_ref = 24.2
# parameter values
Θ1 = 6.28
Θ2 = 25.5
Θ3 = 273
Θ4 = 1.79
Θ5 = 1.75
Θ6 = 1.11
Θ7 = 0.191
Θ8 = 42.3
Θ9 = 9.06
Θ10 = -0.0156
Θ11 = -0.00286
Θ12 = 33.6
Θ13 = -0.0138
Θ14 = 68.3
Θ15 = 2.10
Θ16 = 1.30
Θ17 = 1.42
Θ18 = 0.68
η1 = 1 # change values
η2 = 1
η3 = 1
η4 = 1
η5 = 1
η6 = 1
# f functions (everything ref is related to male, 35 yo, 170 cm, 70 kg)
f_ageing(x) = exp(x*AGE-AGE_ref)
f_ageing_ref(x) = exp(x*AGE_ref-AGE_ref)
f_sigmoid(x, E50, λ) = (x^λ)/(x^λ+E50^λ)
f_central(x) = f_sigmoid(x, Θ12, 1)
f_CLMaturation = f_sigmoid(PMA, Θ8, Θ9)
f_CLMaturation_ref = f_sigmoid(PMA, Θ8, Θ9)
f_Q3Maturation = f_sigmoid(AGE + weeks_40, Θ14, 1)
f_Q3Maturation_ref = f_sigmoid(AGE_ref + weeks_40, Θ14, 1)
f_opiates(x) = 1 # due to absence of f_opiates
f_Al = (1.11 + (1-1.11)/(1 + (AGE/7.1)^(-1.1)))*(9270*WGT)/(8780 + 244*BMI)
f_Al_ref = (0.88 + (1-0.88)/(1 + (AGE/13.4)^(-12.7)))*(9270*WGT)/(6680 + 216*BMI)
# compartments
V1_arterial = Θ1*f_central(WGT)/f_central(WGT_ref)*exp(η1) # [l]
V1_venous = V1_arterial*(1+Θ17*(1-f_central(WGT))) # [l]
V2 = Θ2*WGT/WGT_ref*f_ageing(Θ10)*exp(η2) # [l]
V2_ref = Θ2*WGT_ref/WGT_ref*f_ageing(Θ10)*exp(η2) # [l]
V3 = Θ3*f_Al/f_Al_ref*f_opiates(Θ13)*exp(η3) # [l]
V3_ref = Θ3*f_Al_ref/f_Al_ref*f_opiates(Θ13)*exp(η3) # [l]
CL = Θ15*(WGT/WGT_ref)^(0.75)*(f_CLMaturation/f_CLMaturation_ref)*f_opiates(x)*Θ11*exp(η4) # [l/min], elimination clearance
Q2_arterial = Θ5*(V2/V2_ref)^(0.75)*(1+Θ16*(1-f_Q3Maturation))*exp(η5)
Q2_venous = Q2_arterial*Θ18 # [l/min], compartment clearance
Q3 = Θ6*(V3/V3_ref)^(0.75)*(f_Q3Maturation/f_Q3Maturation_ref)*exp(η6) # [l/min], compartment clearance
# drug transfer rate constants
k10 = CL/V1_venous # [1/s]
k12 = Q2_venous/V1_venous # [1/s]
k13 = Q3/V1_venous # [1/s]
k21 = Q2_venous/V2 # [1/s]
k31 = Q3/V3 # [1/s]
function func(du, u, p, t)
x1, x2, x3, y = u
k10, k12, k13, k21, k31, V1_venous = p
du[1] = -(k12+k12+k13)*x1 + k12*x2 + k13*x3 + (1/V1_venous)*y # dx1
du[2] = k12*x1 - k21*x2 # dx2
du[3] = k31*x1 - k31*x3 # dx3
du[4] = 0 # dy, could be wrong if bolus dose
end
# parameters
p = [k10, k12, k13, k21, k31, V1_venous]
# initial values
u0 = [0.5; 0.25; 0.1; 0.25; 0.25; 1.0]
# timespan
tspan = (0.0,100.0)
# Creating ODE problem and solving the differential equations
prob1 = ODEProblem(func, u0, tspan, p)
sol1 = DifferentialEquations.solve(prob1, saveat=0.1)
plot(sol1)
# Simulating data
odedata = Array(sol1) + 0.01 * randn(size(Array(sol1)))
plot(sol1, alpha = 0.3, legend = false); scatter!(sol1.t, odedata')
# Direct Handling of Bayesian Estimation with Turing
@model function funcModel(data, prob1)
σ = 0.01
# shoud change to uniform distribution
k10 ~ truncated(Normal(0,0.1),0,1)
k12 ~ truncated(Normal(0.22,0.1),0,1)
k13 ~ truncated(Normal(0.1,0.1),0,1)
k21 ~ truncated(Normal(0,0.1),0,1)
k31 ~ truncated(Normal(0,0.1),0,1)
V1_venous ~ truncated(Normal(25,0.5),0,1)
p = [k10, k12, k13, k21, k31, V1_venous]
prob = remake(prob1, p=p)
predicted = DifferentialEquations.solve(prob)
for i = 1:length(predicted)
data[:,i] ~ MvNormal(predicted[i], σ)
end
end
model = funcModel(odedata, prob1)
chain = sample(model, NUTS(.65), MCMCThreads(), 300, 2, progress=true) # crashes big time!
@df chain corrplot([:k10 :k12 :k13 :k2 :k31 :V1_venous], grid = true, bins=50)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment