b1AR system simplification - NRVM CaMKII model

Fitting sensitivity of b1AR system to ISO.

using CurveFit
using DiffEqCallbacks
using Model
using Model: Hz, hil, hilr, second, μM
using ModelingToolkit
using OrdinaryDiffEq
using Plots
using SteadyStateDiffEq
Plots.default(lw=1.5)

Setup b1AR system¶

@parameters ATP = 5000μM ISO = 0μM
@time "Build system" sys = Model.get_bar_sys(ATP, ISO) |> mtkcompile
@time "Build problem" prob = SteadyStateProblem(sys, [])

alg = DynamicSS(KenCarp47())
@time "Solve problem" sol = solve(prob, alg; abstol=1e-10, reltol=1e-10)

Build system: 13.836680 seconds (16.68 M allocations: 1.115 GiB, 4.39% gc time, 99.71% compilation time: <1% of which was recompilation)
Build problem: 5.840855 seconds (7.93 M allocations: 572.274 MiB, 2.00% gc time, 98.80% compilation time: 9% of which was recompilation)
Solve problem: 5.514884 seconds (10.05 M allocations: 693.167 MiB, 3.07% gc time, 99.46% compilation time)

retcode: Success
u: 27-element Vector{Float64}:
 0.02772130832754301
 0.010979484922291234
 0.006296120927385307
 0.005516463242490487
 4.450515338841977
 5.434166833821879
 8.130178795563936
 0.05167662482561576
 0.06164282943482401
 0.036126491312718945
 ⋮
 0.002234274933119501
 0.0014301634098038384
 0.0005920372298763777
 2.2725356124804828e-7
 0.0006214300960438064
 0.008290582144886618
 0.00048603729186110585
 2.291306749472233e-11
 4.679551657897097e-13

@unpack cAMP, fracPKACI, fracPKACII, fracPP1, fracPLBp, fracPLMp, TnI_PKAp, LCCa_PKAp, LCCb_PKAp, IKUR_PKAp, RyR_PKAp = sys
tracked_states = [cAMP, fracPKACI, fracPKACII, fracPP1, fracPLBp, fracPLMp, TnI_PKAp, LCCa_PKAp, LCCb_PKAp, IKUR_PKAp, RyR_PKAp]
stimap = Dict(k=>i for (i, k) in enumerate(tracked_states))

@time "Extract tracked states" sol[tracked_states]

Extract tracked states: 0.634798 seconds (817.57 k allocations: 57.718 MiB, 3.53% gc time, 99.42% compilation time)

11-element Vector{Float64}:
 0.2904076028584033
 0.07298501047186572
 0.1836902046771141
 0.9419363766004317
 0.07669979995815034
 0.11321180903795582
 0.06357879055488538
 0.22065852969961947
 0.2518448370954123
 0.43917939689164937
 0.20534302464846674

Ensemble simulations¶

Log scale for ISO concentration.

iso = exp10.(range(log10(1e-4μM), log10(1μM), length=501))

@time "Solve problems" sim = map(iso) do i
    newprob = remake(prob, p=[ISO => i])
    solve(newprob, alg; abstol=1e-10, reltol=1e-10)[tracked_states]
end

yymat = reduce(hcat, sim)'

Solve problems: 34.887256 seconds (111.39 M allocations: 6.541 GiB, 2.04% gc time, 21.79% compilation time)

501×11 adjoint(::Matrix{Float64}) with eltype Float64:
 0.295162  0.0744678  0.187064  0.941171  …  0.255352  0.443719  0.2085
 0.295249  0.0744951  0.187126  0.941157     0.255416  0.443802  0.208557
 0.295339  0.074523   0.187189  0.941142     0.255481  0.443887  0.208616
 0.29543   0.0745514  0.187253  0.941128     0.255548  0.443972  0.208676
 0.295522  0.0745802  0.187319  0.941113     0.255616  0.44406   0.208737
 0.295617  0.0746096  0.187386  0.941098  …  0.255685  0.444149  0.2088
 0.295713  0.0746396  0.187454  0.941083     0.255755  0.444239  0.208863
 0.295811  0.0746701  0.187523  0.941067     0.255826  0.444332  0.208927
 0.29591   0.0747011  0.187593  0.941051     0.255899  0.444425  0.208993
 0.296012  0.0747327  0.187665  0.941035     0.255973  0.444521  0.20906
 ⋮                                        ⋱                      ⋮
 1.13119   0.270035   0.524632  0.893022     0.495813  0.693572  0.443326
 1.13155   0.270094   0.524707  0.893015     0.495849  0.693603  0.443364
 1.1319    0.270152   0.52478   0.893008     0.495885  0.693633  0.443401
 1.13223   0.270209   0.52485   0.893002  …  0.495919  0.693661  0.443437
 1.13257   0.270265   0.524921  0.892996     0.495953  0.69369   0.443473
 1.1329    0.270321   0.524991  0.892989     0.495987  0.693719  0.443509
 1.13322   0.270374   0.525058  0.892983     0.49602   0.693746  0.443543
 1.13354   0.270428   0.525126  0.892977     0.496053  0.693774  0.443577
 1.13385   0.27048    0.525191  0.892971  …  0.496085  0.6938    0.44361

xopts = (xlims=(iso[begin], iso[end]), minorgrid=true, xscale=:log10, xlabel="ISO (μM)",)
plot(iso, yymat[:, stimap[cAMP]]; lab="cAMP", ylabel="Conc. (μM)", legend=:topleft, xopts...)

plot(iso, yymat[:, stimap[fracPKACI]]; lab="PKACI", ylabel="Activation fraction")
plot!(iso, yymat[:, stimap[fracPKACII]], lab="PKACII")
plot!(iso, yymat[:, stimap[fracPP1]], lab="PP1", legend=:topleft; xopts...)

PKACI activity¶

pkaci_model(p, x) = @. p[1] * hil(x, p[2]) + p[3]
xdata = iso
ydata = yymat[:, stimap[fracPKACI]]
p0 = [0.3, 0.01μM, 0.08]
prob = NonlinearCurveFitProblem(pkaci_model, p0, xdata, ydata)
@time "Fit PKACI" fit_pkac1 = solve(prob)

println("PKACI")
println("Basal activity: ", fit_pkac1.u[3])
println("Activated activity: ", fit_pkac1.u[1])
println("Michaelis constant: ", fit_pkac1.u[2], " μM")
println("RMSE: ", mse(fit_pkac1) |> sqrt)

plot(xdata, [ydata predict(fit_pkac1)], lab=["Full model" "Fitted"], line=[:dash :dot], title="PKACI", legend=:topleft; xopts...)

Fit PKACI: 4.491757 seconds (4.59 M allocations: 320.706 MiB, 1.48% gc time, 99.93% compilation time)
PKACI
Basal activity: 0.0734020796431307
Activated activity: 0.1993769973447233
Michaelis constant: 0.013906485023051387 μM
RMSE: 0.00032795425080978406

savefig("pkaci_fit.pdf")

"/home/runner/work/camkii-cardiomyocyte-model/camkii-cardiomyocyte-model/docs/fits/pkaci_fit.pdf"

plot(xdata, residuals(fit_pkac1) ./ ydata .* 100; title="PKACI error (%)", lab=false, xopts...)

PKACII activity¶

pkacii_model(p, x) = @. p[1] * hil(x, p[2]) + p[3]
xdata = iso
ydata = yymat[:, stimap[fracPKACII]]
p0 = [0.4, 0.01μM, 0.2]
prob = NonlinearCurveFitProblem(pkacii_model, p0, xdata, ydata)
@time fit_pkac2 = solve(prob)

println("PKACII")
println("Basal activity: ", fit_pkac2.u[3])
println("Activated activity: ", fit_pkac2.u[1])
println("Michaelis constant: ", fit_pkac2.u[2], " μM")
println("RMSE: ", mse(fit_pkac2) |> sqrt)

plot(xdata, [ydata predict(fit_pkac2)], lab=["Full model" "Fitted"], line=[:dash :dot], title="PKACII", legend=:topleft; xopts...)

  3.528428 seconds (3.65 M allocations: 256.746 MiB, 2.01% gc time, 99.93% compilation time)
PKACII
Basal activity: 0.1839918426023136
Activated activity: 0.34437525651872514
Michaelis constant: 0.01025117760817389 μM
RMSE: 0.00024441507138846913

savefig("pkacii_fit.pdf")

"/home/runner/work/camkii-cardiomyocyte-model/camkii-cardiomyocyte-model/docs/fits/pkacii_fit.pdf"

plot(xdata, residuals(fit_pkac2) ./ ydata .* 100; title="PKACII error (%)", lab=false, xopts...)

PP1 activity¶

pp1_model(p, x) = @. p[1] * hil(p[2], x) + p[3]
xdata = iso
ydata = yymat[:, stimap[fracPP1]]
p0 = [0.1, 3e-3μM, 0.8]
prob = NonlinearCurveFitProblem(pp1_model, p0, xdata, ydata)
@time fit_pp1 = solve(prob)

println("PP1")
println("Repressible activity: ", fit_pp1.u[1])
println("Minimal activity: ", fit_pp1.u[3])
println("Repressive Michaelis constant: ", fit_pp1.u[2], " μM")
println("RMSE: ", mse(fit_pp1) |> sqrt)

plot(xdata, [ydata predict(fit_pp1)], lab=["Full model" "Fitted"], line=[:dash :dot], title="PP1", legend=:topright; xopts...)

  3.541570 seconds (3.78 M allocations: 266.081 MiB, 1.36% gc time, 99.93% compilation time)
PP1
Repressible activity: 0.04919243277385461
Minimal activity: 0.8927127483276253
Repressive Michaelis constant: 0.006361978464789555 μM
RMSE: 3.551130082206751e-5

savefig("pp1_fit.pdf")

"/home/runner/work/camkii-cardiomyocyte-model/camkii-cardiomyocyte-model/docs/fits/pp1_fit.pdf"

plot(xdata, residuals(fit_pp1) ./ ydata .* 100; title="PP1 error (%)", lab=false, xopts...)

PLBp¶

xdata = iso
ydata = yymat[:, stimap[fracPLBp]]
plot(xdata, ydata, title="PLBp fraction", lab=false; xopts...)

plb_model(p, x) = @. p[1] * hil(x, p[2], p[3]) + p[4]
p0 = [0.8, 1e-2μM, 1.0, 0.1]
prob = NonlinearCurveFitProblem(plb_model, p0, xdata, ydata)
@time fit_plb = solve(prob)

println("PLBp")
println("Basal activity: ", fit_plb.u[4])
println("Activated activity: ", fit_plb.u[1])
println("Michaelis constant: ", fit_plb.u[2], " μM")
println("Hill coefficient: ", fit_plb.u[3])
println("RMSE: ", mse(fit_plb) |> sqrt)

plot(xdata, [ydata predict(fit_plb)], lab=["Full model" "Fitted"], line=[:dash :dot], title="PLBp", legend=:topleft; xopts...)

  3.613462 seconds (3.90 M allocations: 274.996 MiB, 2.24% gc time, 99.67% compilation time)
PLBp
Basal activity: 0.08300162154344362
Activated activity: 0.7923506041631098
Michaelis constant: 0.005936779858771097 μM
Hill coefficient: 1.8385090667075448
RMSE: 0.008452352383475832

savefig("plbp_fit.pdf")

"/home/runner/work/camkii-cardiomyocyte-model/camkii-cardiomyocyte-model/docs/fits/plbp_fit.pdf"

plot(xdata, residuals(fit_plb) ./ ydata .* 100; title="PLBp error (%)", lab=false, xopts...)

PLMp¶

xdata = iso
ydata = yymat[:, stimap[fracPLMp]]
plot(xdata, ydata, title="PLMp fraction", lab=false; xopts...)

plm_model(p, x) = @. p[1] * hil(x, p[2], p[3]) + p[4]
p0 = [0.8, 1e-2μM, 1.0, 0.1]
prob = NonlinearCurveFitProblem(plm_model, p0, xdata, ydata)
@time fit_plm = solve(prob)

println("PLMp")
println("Basal activity: ", fit_plm.u[4])
println("Activated activity: ", fit_plm.u[1])
println("Michaelis constant: ", fit_plm.u[2], " μM")
println("Hill coefficient: ", fit_plm.u[3])
println("RMSE: ", mse(fit_plm) |> sqrt)

plot(xdata, [ydata predict(fit_plm)], lab=["Full model" "Fitted"], line=[:dash :dot], title="PLMp", legend=:topleft; xopts...)

  3.786958 seconds (3.70 M allocations: 261.055 MiB, 5.91% gc time, 99.71% compilation time)
PLMp
Basal activity: 0.11771895584233232
Activated activity: 0.6617431876526159
Michaelis constant: 0.008181899156781835 μM
Hill coefficient: 1.3709432030651023
RMSE: 0.003287961010974935

savefig("plmp_fit.pdf")

"/home/runner/work/camkii-cardiomyocyte-model/camkii-cardiomyocyte-model/docs/fits/plmp_fit.pdf"

plot(xdata, residuals(fit_plm) ./ ydata .* 100; title="PLMp error (%)", lab=false, xopts...)

TnIp¶

xdata = iso
ydata = yymat[:, stimap[TnI_PKAp]]
plot(xdata, ydata, title="TnIp fraction", lab=false; xopts...)

tni_model(p, x) = @. p[1] * hil(x, p[2], p[3]) + p[4]
p0 = [0.8, 1e-2μM, 1.0, 0.1]
lb = [0.1, 1e-9μM, 1.0, 0.0]
prob = NonlinearCurveFitProblem(tni_model, p0, xdata, ydata)
@time fit_tni = solve(prob)

println("TnIp")
println("Basal activity: ", fit_tni.u[4])
println("Activated activity: ", fit_tni.u[1])
println("Michaelis constant: ", fit_tni.u[2], " μM")
println("Hill coefficient: ", fit_tni.u[3])
println("RMSE: ", mse(fit_tni) |> sqrt)

plot(xdata, [ydata predict(fit_tni)], lab=["Full model" "Fitted"], line=[:dash :dot], title="TnIp", legend=:topleft; xopts...)

  3.620519 seconds (3.70 M allocations: 261.413 MiB, 2.20% gc time, 99.68% compilation time)
TnIp
Basal activity: 0.06752245015254726
Activated activity: 0.7482251059841436
Michaelis constant: 0.007856634350569719 μM
Hill coefficient: 1.6971951371664125
RMSE: 0.007467040874206868

savefig("tni_fit.pdf")

"/home/runner/work/camkii-cardiomyocyte-model/camkii-cardiomyocyte-model/docs/fits/tni_fit.pdf"

plot(xdata, residuals(fit_tni) ./ ydata .* 100; title="TnIp error (%)", lab=false, xopts...)

LCCap¶

xdata = iso
ydata = yymat[:, stimap[LCCa_PKAp]]
plot(xdata, ydata, title="LCCap fraction", lab=false; xopts...)

lcc_model(p, x) = @. p[1] * hil(x, p[2]) + p[3]
p0 = [0.8, 1e-2μM, 0.1]
prob = NonlinearCurveFitProblem(lcc_model, p0, xdata, ydata)
@time fit_lcca = solve(prob)

println("LCCap")
println("Basal activity: ", fit_lcca.u[3])
println("Activated activity: ", fit_lcca.u[1])
println("Michaelis constant: ", fit_lcca.u[2], " μM")
println("RMSE: ", mse(fit_lcca) |> sqrt)

plot(xdata, [ydata predict(fit_lcca)], lab=["Full model" "Fitted"], line=[:dash :dot], title="LCCap", legend=:topleft; xopts...)

  3.663281 seconds (3.65 M allocations: 257.043 MiB, 1.71% gc time, 99.92% compilation time)
LCCap
Basal activity: 0.22081895230233178
Activated activity: 0.23338438387766694
Michaelis constant: 0.007263438339117075 μM
RMSE: 0.00013512020520144722

savefig("lcca_fit.pdf")

"/home/runner/work/camkii-cardiomyocyte-model/camkii-cardiomyocyte-model/docs/fits/lcca_fit.pdf"

plot(xdata, residuals(fit_lcca) ./ ydata .* 100; title="LCCap error (%)", lab=false, xopts...)

LCCbp¶

xdata = iso
ydata = yymat[:, stimap[LCCb_PKAp]]
plot(xdata, ydata, title="LCCbp fraction", lab=false; xopts...)

p0 = [0.8, 1e-2μM, 0.1]
lb = [0.1, 1e-9μM, 0.0]
prob = NonlinearCurveFitProblem(lcc_model, p0, xdata, ydata)
@time fit_lccb = solve(prob)

println("LCCbp")
println("Basal activity: ", fit_lccb.u[3])
println("Activated activity: ", fit_lccb.u[1])
println("Michaelis constant: ", fit_lccb.u[2], " μM")
println("RMSE: ", mse(fit_lccb) |> sqrt)

plot(xdata, [ydata predict(fit_lccb)], lab=["Full model" "Fitted"], line=[:dash :dot], title="LCCbp", legend=:topleft; xopts...)

  0.000945 seconds (784 allocations: 2.744 MiB)
LCCbp
Basal activity: 0.25200744727086005
Activated activity: 0.24559143115771226
Michaelis constant: 0.006961207541722576 μM
RMSE: 0.000138464877691979

savefig("lccbp_fit.pdf")

"/home/runner/work/camkii-cardiomyocyte-model/camkii-cardiomyocyte-model/docs/fits/lccbp_fit.pdf"

plot(xdata, residuals(fit_lccb) ./ ydata .* 100; title="LCCbp error (%)", lab=false, xopts...)

KURp¶

xdata = iso
ydata = yymat[:, stimap[IKUR_PKAp]]
plot(xdata, ydata, title="KURp fraction", lab=false; xopts...)

kur_model(p, x) = @. p[1] * hil(x, p[2]) + p[3]
p0 = [0.8, 1e-2μM, 0.1]
lb = [0.1, 1e-9μM, 0.0]
prob = NonlinearCurveFitProblem(kur_model, p0, xdata, ydata)
@time fit_kur = solve(prob)

println("KURp")
println("Basal activity: ", fit_kur.u[3])
println("Activated activity: ", fit_kur.u[1])
println("Michaelis constant: ", fit_kur.u[2], " μM")
println("RMSE: ", mse(fit_kur) |> sqrt)

plot(xdata, [ydata predict(fit_kur)], lab=["Full model" "Fitted"], line=[:dash :dot], title="KURp", legend=:topleft; xopts...)

  3.616514 seconds (3.65 M allocations: 256.945 MiB, 1.66% gc time, 99.92% compilation time)
KURp
Basal activity: 0.43936332364813485
Activated activity: 0.2556595334759436
Michaelis constant: 0.005578209608528763 μM
RMSE: 0.00014894273014033873

savefig("kurp_fit.pdf")

"/home/runner/work/camkii-cardiomyocyte-model/camkii-cardiomyocyte-model/docs/fits/kurp_fit.pdf"

plot(xdata, residuals(fit_kur) ./ ydata .* 100; title="KURp error (%)", lab=false, xopts...)

RyRp¶

xdata = iso
ydata = yymat[:, stimap[RyR_PKAp]]
plot(xdata, ydata, title="RyRp fraction", lab=false; xopts...)

ryr_model(p, x) = @. p[1] * x / (x + p[2]) + p[3]
p0 = [0.3, 1e-2μM, 0.1]
lb = [0.0, 1e-9μM, 0.0]
prob = NonlinearCurveFitProblem(ryr_model, p0, xdata, ydata)
@time fit_ryr = solve(prob)

println("RyRp")
println("Basal activity: ", fit_ryr.u[3])
println("Activated activity: ", fit_ryr.u[1])
println("Michaelis constant: ", fit_ryr.u[2], " μM")
println("RMSE: ", mse(fit_ryr) |> sqrt)

plot(xdata, [ydata predict(fit_ryr)], lab=["Full model" "Fitted"], line=[:dash :dot], title="RyRp", legend=:topleft; xopts...)

  3.682099 seconds (3.87 M allocations: 271.534 MiB, 1.61% gc time, 99.93% compilation time)
RyRp
Basal activity: 0.20539770176393893
Activated activity: 0.23988888631204675
Michaelis constant: 0.007509969338942928 μM
RMSE: 7.705670671540831e-5

savefig("ryrp_fit.pdf")

"/home/runner/work/camkii-cardiomyocyte-model/camkii-cardiomyocyte-model/docs/fits/ryrp_fit.pdf"

plot(xdata, residuals(fit_ryr) ./ ydata .* 100; title="RyRp error (%)", lab=false, xopts...)

This notebook was generated using Literate.jl.