-
Notifications
You must be signed in to change notification settings - Fork 17
/
practical_identifiability.jl
115 lines (100 loc) · 3.71 KB
/
practical_identifiability.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
#import Pkg; Pkg.add("DifferentialEquations")
#import Pkg; Pkg.add("ForwardDiff")
#import Pkg; Pkg.add("DiffResults")
#import Pkg; Pkg.add("DiffEqSensitivity")
#import Pkg; Pkg.add("DataFrames")
#import Pkg; Pkg.add("CSV")
#import Pkg; Pkg.add("Dierckx")
#import Pkg; Pkg.add("StatsPlots")
#import Pkg; Pkg.add("LaTeXStrings")
#import Pkg; Pkg.add("Plots")
using DifferentialEquations
using ForwardDiff
using DiffResults
using DiffEqSensitivity
using Statistics
using Random
using LinearAlgebra
using DataFrames
using CSV
using Dierckx
using Plots
using StatsPlots
using LaTeXStrings
function f_apop(dx, x, p, t)
x1, x2, x3, x4, x5, x6 = x
E, tp, ti, td, k, Rm, a1, C1, C2, C4, C5, Ub, U0, Um, Rg, alpha, beta = p
Vp = 3
Vi = 11
Vg = 10
C3 = 100/100
meal_t = 300, 650, 1100
meal_q = 60e3, 40e3, 50e3
f1 = Rm / (1 + exp(x3 / (Vg * C1) - a1))
f2 = Ub * (1 - exp(-x3 / (Vg * C2)))
kappa = (1 / Vi + 1 / (E * ti)) / C4
f3 = (U0 + Um / (1 + (kappa * x2)^(-beta))) / (Vg * C3)
f4 = Rg / (1 + exp(alpha * (1 - x6 / (Vp * C5))))
dt1 = t - meal_t[1]
dt2 = t - meal_t[2]
dt3 = t - meal_t[3]
IG1 = 0.5 * meal_q[1] * k * exp(-k * dt1) * (sign(dt1) + 1)
IG2 = 0.5 * meal_q[2] * k * exp(-k * dt2) * (sign(dt2) + 1)
IG3 = 0.5 * meal_q[3] * k * exp(-k * dt3) * (sign(dt3) + 1)
IG = IG1 + IG2 + IG3
tmp = E * (x1 / Vp - x2 / Vi)
dx[1] = (f1 - tmp - x1 / tp)
dx[2] = (tmp - x2 / ti)
dx[3] = (f4 + IG - f2 - f3 * x3)
dx[4] = (x1 - x4) / td
dx[5] = (x4 - x5) / td
dx[6] = (x5 - x6) / td
end
cscale = 1; tscale = 1.;
p = [0.200853979,5.986360211,101.2036365,11.97776077,0.00833154,208.6221286,6.592132806,301.2623884,37.65196466,78.75855574,25.93801618,71.32611075,4.063671595,88.97766871,179.85672,7.536059566,1.783340558]
x0 = [36., 44., 11000., 0., 0., 0.] / cscale
tspan = (0.0, 1800.0)
prob_apop = ODELocalSensitivityProblem(f_apop, x0, tspan, p)
sol_apop = solve(prob_apop, alg_hints=[:stiff], saveat=0.1)
x_apop, dp_apop = extract_local_sensitivities(sol_apop)
lab = [L"E", L"tp", L"ti", L"td", L"k", L"Rm", L"a1", L"C1", L"C2", L"C4", L"C5", L"Ub", L"U0", L"Um", L"Rg", L"alpha", L"beta"]
σ = 0.01 * std(x_apop, dims=2)
cov_ϵ = σ[4]
dp = dp_apop
cols = 4:4
#plot(sol_apop.t, x_apop[4,:], lw=2)
Nt = length(dp[1][1,:])
Nstate = length(dp[1][:,1])
Nparam = length(dp[:,1])
F = zeros(Float64, Nparam, Nparam)
perm = vcat(1, sort(rand(2:Nt-1, Nt÷5)), Nt)
for i in perm
S = reshape(dp[1][:,i], (Nstate,1))
for j = 2:Nparam
S = hcat(S, reshape(dp[j][:,i], (Nstate,1)))
end
global F += S[cols,:]' * inv(cov_ϵ) * S[cols,:]
end
C = inv(F)
R = ones(size(C))
R = [C[i,j]/sqrt(C[i,i]*C[j,j]) for i = 1:size(C)[1], j = 1:size(C)[1]]
#heatmap(R, xlims=(0.5,size(R)[1]+0.5), aspect_ratio = 1, color = :inferno, clims = (-1, 1),
# xticks = (1:1:size(C)[1], lab), xtickfont = font(14, "Times"),
# yticks = (1:1:size(C)[1], lab), ytickfont = font(14, "Times"), fmt = :pdf, dpi=300)
#savefig("correlation_matrix")
abs.(R) .> 0.99
lowerbound = sqrt.(diag(inv(F))) / tscale
lowerbound[1:3:7] = lowerbound[1:3:7] / cscale
for i = 1:length(lab)
println(lab[i], '\t', lowerbound[i])
end
for i = 1:Nparam
println(eigvals(F)[i])
println(eigvecs(F)[:,i])
println('\n')
end
#bar(eigvecs(F)[:,1:9], ylabel = "FIM null eigenvector coefficients", ytickfont = font(12, "Times"),
# xticks = (1:1:size(C)[1], lab), xtickfont = font(12, "Courier"),
# legendfontsize = 10, label = [L"null_1" L"null_2" L"null_3" L"null_4" L"null_5" L"null_6" L"null_7" L"null_8" L"null_9"], fmt = :png,
# legend=:topright, dpi=300)
#savefig("nulleigen_apop")