using ModelingToolkit
using OrdinaryDiffEq, SteadyStateDiffEq, DiffEqCallbacks
using Plots
using CurveFit
using Model
using Model: μM, hil, second, Hz, get_camkii_sys
Plots.default(lw=1.5)CaMKII sensitivity to calcium¶
Modeling CaM-Cax binding to CaMKII only. No phosphorylation or oxidation reactions.
@parameters Ca = 0μM ROS = 0μM
@time "Build system" sys = Model.get_camkii_sys(; Ca=Ca, ROS=ROS) |> mtkcompile
@time "Build problem" prob = SteadyStateProblem(sys, [sys.k_phosCaM => 0])Build system: 9.931878 seconds (12.18 M allocations: 847.185 MiB, 2.29% gc time, 99.24% compilation time: <1% of which was recompilation)
Build problem: 7.368647 seconds (11.79 M allocations: 846.155 MiB, 2.73% gc time, 98.86% compilation time: 8% of which was recompilation)
SteadyStateProblem with uType Vector{Float64}. In-place: true
u0: 17-element Vector{Float64}:
0.0
0.0
0.0
3.3
0.0
0.0
0.01581
0.010750000000000001
1.01
0.62372
0.00666
0.0166
0.39868000000000003
1.0
2.0e-5
0.008577999999999999
0.08433Physiological cytosolic calcium levels ranges from 30nM to 10μM.
ca = exp10.(range(log10(0.03μM), log10(10μM), length=1001))
prob_func = (prob, i, repeat) -> (prob.ps[Ca] = ca[i]; prob)
alg = DynamicSS(TRBDF2())
sol0 = solve(prob, alg) ## warmup
@time "Solve problem" sim = solve(EnsembleProblem(prob; prob_func), alg; trajectories=length(ca), abstol=1e-8, reltol=1e-8);Solve problem: 16.320868 seconds (41.98 M allocations: 9.994 GiB, 16.38% gc time, 35.83% compilation time)
"""Extract values from ensemble simulations by a symbol"""
extract(sim, k) = map(s -> s[k], sim)Main.var"##225".extractStatus of the CaMKII system across a range of calcium concentrations.
xopts = (xlabel="Ca (μM)", xscale=:log10, minorgrid=true, xlims=(ca[1], ca[end]))
let
plot(ca, extract(sim, sys.CaM0), lab="CaM0", ylabel="Conc. (μM)"; xopts...)
plot!(ca, extract(sim, sys.Ca2CaM_C), lab="Ca2CaM_C")
plot!(ca, extract(sim, sys.Ca2CaM_N), lab="Ca2CaM_N")
plot!(ca, extract(sim, sys.Ca4CaM), lab="Ca4CaM")
plot!(ca, extract(sim, sys.CaMK), lab="CaMK")
plot!(ca, extract(sim, sys.CaM0_CaMK), lab="CaM0_CaMK")
plot!(ca, extract(sim, sys.Ca2CaM_C_CaMK), lab="Ca2CaM_C_CaMK")
plot!(ca, extract(sim, sys.Ca2CaM_N_CaMK), lab="Ca2CaM_N_CaMK")
plot!(ca, extract(sim, sys.Ca4CaM_CaMK), lab="Ca4CaM_CaMK", legend=:topleft)
end
savefig("camkii_cam_binding.pdf")"/home/runner/work/camkii-cardiomyocyte-model/camkii-cardiomyocyte-model/docs/fits/camkii_cam_binding.pdf"CaMKAct = 1 - (sys.CaMK + sys.CaM0_CaMK) / sys.CAMKII_T
println("Basal activity with 30nM Ca is ", sim[1][CaMKAct])
xdata = ca
ydata = extract(sim, CaMKAct)
plot(xdata, ydata, label=false, title="Active CaMKII", ylims=(0, 0.5); xopts...)Basal activity with 30nM Ca is 0.0005733143444885958
Least-square fitting¶
Least-square fitting of steady state CaMKII activities against calcium concentrations.
Old model¶
Using a single Hill function to fit steady state CaMKII activities.
model_camk_single_hill(p, x) = @. p[1] * hil(x, p[2], p[4]) + p[3]
p0 = [0.4, 1μM, 0.0, 2.0]
prob = NonlinearCurveFitProblem(model_camk_single_hill, p0, xdata, ydata)
@time fit = solve(prob)
println("Basal activity: ", fit.u[3])
println("Maximal activated activity: ", fit.u[1])
println("Half-saturation Ca concentration: ", fit.u[2], " μM")
println("Hill coefficient: ", fit.u[4])
println("RMSE: ", mse(fit) |> sqrt) 4.873181 seconds (4.73 M allocations: 333.641 MiB, 99.58% compilation time)
Basal activity: 0.0051912756281683485
Maximal activated activity: 0.42996709609968464
Half-saturation Ca concentration: 0.9630642784842367 μM
Hill coefficient: 2.298409070803075
RMSE: 0.004835754512986751
plot(xdata, [ydata predict(fit)], lab=["Full model" "Fitted"], line=[:dash :dot], title="Single Hill function fit", legend=:topleft, ylabel="Bound CaMKII"; xopts...)
model_camk(p, x) = @. p[1] * hil(x, p[2], 2) + p[3] * hil(x, p[4], 4) + p[5]
p0 = [0.4, 1μM, 0.2, 1μM, 0.0]
prob = NonlinearCurveFitProblem(model_camk, p0, xdata, ydata)
@time fit = solve(prob)
println("Basal activity: ", fit.u[5])
println("Maximal activity by CaM-Ca2 binding: ", fit.u[1])
println("Half saturation Ca concentration for CaM-Ca2 binding: ", fit.u[2], " μM")
println("Maximal activity by CaM-Ca4 binding: ", fit.u[3])
println("Half saturation Ca concentration for CaM-Ca4 binding: ", fit.u[4], " μM")
println("RMSE: ", mse(fit) |> sqrt) 4.070747 seconds (4.42 M allocations: 315.271 MiB, 99.77% compilation time)
Basal activity: 0.0010196222878298715
Maximal activity by CaM-Ca2 binding: 0.26511806820971623
Half saturation Ca concentration for CaM-Ca2 binding: 0.738525927129968 μM
Maximal activity by CaM-Ca4 binding: 0.1635740292981544
Half saturation Ca concentration for CaM-Ca4 binding: 1.2514941370204529 μM
RMSE: 0.0009703936611898563
plot(xdata, [ydata predict(fit)], lab=["Full model" "Fitted"], line=[:dash :dot], title="Dual Hill function fit", legend=:topleft, ylabel="Bound CaMKII"; xopts...)
savefig("camkii_act_fit.pdf")"/home/runner/work/camkii-cardiomyocyte-model/camkii-cardiomyocyte-model/docs/fits/camkii_act_fit.pdf"Polynomial fitting¶
prob = CurveFitProblem(xdata, ydata)
alg = RationalPolynomialFitAlgorithm(num_degree=4, den_degree=4)
@time sol = solve(prob, alg)
println("Numerator coefficients: ", sol.u[1:5])
println("Denominator coefficients: ", vcat(1.0, sol.u[6:end]))
println("RMSE: ", mse(sol) |> sqrt) 4.884846 seconds (5.60 M allocations: 394.114 MiB, 99.79% compilation time)
Numerator coefficients: [0.0004422630332169665, -0.019898724925566106, 0.9341965745698904, -0.1638740437349248, 0.7550657908815218]
Denominator coefficients: [1.0, 2.9612109541798524, 1.2704507719309002, -0.19362315594584553, 1.7553801446919701]
RMSE: 4.393934501487125e-5
plot(xdata, [ydata sol.(xdata)], lab=["Full model" "Fitted"], line=[:dash :dot], title="Rational polynomial fit", legend=:topleft, ylabel="Bound CaMKII"; xopts...)
This notebook was generated using Literate.jl.