Correct camkii params (#37)

* pkg

* rm 3Hz

* internalize phos rate

* wip

* pkg

* wip

Manifest.toml

@@ -382,9 +382,9 @@ version = "1.9.1"
 
 [[deps.DiffEqBase]]
 deps = ["ArrayInterface", "ConcreteStructs", "DataStructures", "DocStringExtensions", "EnumX", "EnzymeCore", "FastBroadcast", "FastClosures", "FastPower", "ForwardDiff", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "Logging", "Markdown", "MuladdMacro", "Parameters", "PreallocationTools", "PrecompileTools", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SciMLStructures", "Setfield", "Static", "StaticArraysCore", "Statistics", "TruncatedStacktraces"]
-git-tree-sha1 = "697abdf4af0e38199e9eabff6ccdf65255de855d"
+git-tree-sha1 = "b7dbeaa770bad0980ddddf606de814cff2acb3bc"
 uuid = "2b5f629d-d688-5b77-993f-72d75c75574e"
-version = "6.159.0"
+version = "6.160.0"
 
     [deps.DiffEqBase.extensions]
     DiffEqBaseCUDAExt = "CUDA"
@@ -1052,9 +1052,9 @@ version = "0.1.17"
 
 [[deps.LazyArrays]]
 deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra", "MacroTools", "SparseArrays"]
-git-tree-sha1 = "360f6039babd6e4d6364eff0d4fc9120834a2d9a"
+git-tree-sha1 = "376bc148ae72e68a08f0d5d8a69e287025a37687"
 uuid = "5078a376-72f3-5289-bfd5-ec5146d43c02"
-version = "2.2.1"
+version = "2.2.2"
 
     [deps.LazyArrays.extensions]
     LazyArraysBandedMatricesExt = "BandedMatrices"
@@ -1404,9 +1404,9 @@ version = "4.1.0"
 
 [[deps.NonlinearSolveBase]]
 deps = ["ADTypes", "Adapt", "ArrayInterface", "CommonSolve", "Compat", "ConcreteStructs", "DifferentiationInterface", "EnzymeCore", "FastClosures", "FunctionProperties", "LinearAlgebra", "Markdown", "MaybeInplace", "Preferences", "Printf", "RecursiveArrayTools", "SciMLBase", "SciMLJacobianOperators", "SciMLOperators", "StaticArraysCore", "SymbolicIndexingInterface", "TimerOutputs"]
-git-tree-sha1 = "545e3a49c20b453a929951657f039716216cdf3c"
+git-tree-sha1 = "46772fc296d9f16c3ab78a8ef00008ab075de677"
 uuid = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
-version = "1.3.2"
+version = "1.3.3"
 
     [deps.NonlinearSolveBase.extensions]
     NonlinearSolveBaseBandedMatricesExt = "BandedMatrices"
@@ -1530,9 +1530,9 @@ version = "1.1.2"
 
 [[deps.OrdinaryDiffEqCore]]
 deps = ["ADTypes", "Accessors", "Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "DocStringExtensions", "EnumX", "FastBroadcast", "FastClosures", "FastPower", "FillArrays", "FunctionWrappersWrappers", "InteractiveUtils", "LinearAlgebra", "Logging", "MacroTools", "MuladdMacro", "Polyester", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SciMLStructures", "SimpleUnPack", "Static", "StaticArrayInterface", "StaticArraysCore", "SymbolicIndexingInterface", "TruncatedStacktraces"]
-git-tree-sha1 = "be2c628185fcf94b544fba86b9722be40c3ac305"
+git-tree-sha1 = "6f86041d5978ce3574f022f4098e0c216e2a3395"
 uuid = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
-version = "1.10.1"
+version = "1.10.2"
 weakdeps = ["EnzymeCore"]
 
     [deps.OrdinaryDiffEqCore.extensions]
@@ -2336,9 +2336,9 @@ version = "7.7.0+0"
 
 [[deps.SymbolicIndexingInterface]]
 deps = ["Accessors", "ArrayInterface", "RuntimeGeneratedFunctions", "StaticArraysCore"]
-git-tree-sha1 = "20cf607cafb31f922bce84d60379203e7a126911"
+git-tree-sha1 = "6c6761e08bf5a270905cdd065be633abfa1b155b"
 uuid = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
-version = "0.3.34"
+version = "0.3.35"
 
 [[deps.SymbolicLimits]]
 deps = ["SymbolicUtils"]
@@ -2362,9 +2362,9 @@ version = "3.7.2"
 
 [[deps.Symbolics]]
 deps = ["ADTypes", "ArrayInterface", "Bijections", "CommonWorldInvalidations", "ConstructionBase", "DataStructures", "DiffRules", "Distributions", "DocStringExtensions", "DomainSets", "DynamicPolynomials", "IfElse", "LaTeXStrings", "Latexify", "Libdl", "LinearAlgebra", "LogExpFunctions", "MacroTools", "Markdown", "NaNMath", "PrecompileTools", "Primes", "RecipesBase", "Reexport", "RuntimeGeneratedFunctions", "SciMLBase", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArraysCore", "SymbolicIndexingInterface", "SymbolicLimits", "SymbolicUtils", "TermInterface"]
-git-tree-sha1 = "0caef7687abf7094132fa3112bf5514c36a99226"
+git-tree-sha1 = "24e006074ef13894ed23d006f55e6082998c9035"
 uuid = "0c5d862f-8b57-4792-8d23-62f2024744c7"
-version = "6.18.3"
+version = "6.19.0"
 
     [deps.Symbolics.extensions]
     SymbolicsForwardDiffExt = "ForwardDiff"

docs/bcl.jl

@@ -49,31 +49,13 @@ plot(sol2, idxs=[sys.Cai_sub_SR, sys.Cai_sub_SL, sys.Cai_mean], tspan=(299, 300)
 #---
 plot(sol2, idxs=sys.CaMKAct, title="Active CaMKII")
 
-# ## 3Hz
-callback = build_stim_callbacks(Istim, stimend; period=1/3, starttime=stimstart)
-sol3 = solve(prob, alg; callback, abstol=1e-6, reltol=1e-6)
-
-#---
-plot(sol3, idxs=sys.vm*1000, title="Action potential")
-
-#---
-plot(sol3, idxs=sys.vm*1000, title="Action potential", tspan=(299, 300))
-
-#---
-plot(sol3, idxs=[sys.Cai_sub_SR, sys.Cai_sub_SL, sys.Cai_mean], tspan=(299, 300), title="Calcium transient")
-
-#---
-plot(sol3, idxs=sys.CaMKAct, title="Active CaMKII")
-
-# ## Comparing 1-3 Hz
-plot(sol, idxs=sys.vm*1000, title="Action potential", lab="1Hz")
-plot!(sol2, idxs=sys.vm*1000, lab="2Hz")
-plot!(sol3, idxs=sys.vm*1000, lab="3Hz", tspan=(299, 300), xlabel="Time (sec.)", ylabel="Voltage (mV)")
+# ## Comparing 1-2 Hz
+plot(sol, idxs=sys.vm*1000, title="Action potential", lab="1Hzf")
+plot!(sol2, idxs=sys.vm*1000, lab="2Hz", tspan=(299, 300), xlabel="Time (sec.)", ylabel="Voltage (mV)")
 
 #---
 plot(sol, idxs=sys.CaMKAct, title="CaMKII", lab="1Hz")
-plot!(sol2, idxs=sys.CaMKAct, lab="2Hz")
-plot!(sol3, idxs=sys.CaMKAct, lab="3Hz", ylim=(0.0, 1.0), xlabel="Time (sec.)", ylabel="Act. fraction (AU)")
+plot!(sol2, idxs=sys.CaMKAct, lab="2Hz",  ylim=(0.0, 1.0), xlabel="Time (sec.)", ylabel="Act. fraction (AU)")
 
 # ## Data fitting
 ### Pacing duration and CaMKII activity

docs/camk-sens.jl

@@ -12,18 +12,19 @@ Plots.default(lw=1.5)
 # ## Setup CaMKII system
 # Model CaM/Calcium binding only. NO phosphorylation or oxidation.
 @parameters Ca = 0μM ROS = 0μM
-sys = get_camkii_sys(Ca; ROS, phospho_rate=0, simplify=true)
+sys = get_camkii_sys(Ca; ROS, simplify=true)
 equations(sys)
 
 #---
-prob = SteadyStateProblem(sys, [])
+prob = SteadyStateProblem(sys, [sys.k_phosCaM => 0])
 alg = DynamicSS(Rodas5P())
-ca = exp10.(range(log10(0.1μM), log10(100μM), length=101))
+logrange(0.01μM, 10μM, length=101)
+ca = logrange(0.03μM, 10μM, length=101)
 prob_func = (prob, i, repeat) -> begin
     remake(prob, p=[Ca => ca[i]])
 end
 trajectories = length(ca)
-sol0 = solve(remake(prob, p=[Ca => 0.0]), alg) ## warmup
+sol0 = solve(prob, alg) ## warmup
 sim = solve(EnsembleProblem(prob; prob_func, safetycopy=false), alg; trajectories)
 
 #---
@@ -42,19 +43,19 @@ plot!(ca .* 1000, extract(sim, sys.CaMK), lab="CaMK")
 plot!(ca .* 1000, extract(sim, sys.CaM0_CaMK), lab="CaM0_CaMK")
 plot!(ca .* 1000, extract(sim, sys.Ca2CaM_C_CaMK), lab="Ca2CaM_C_CaMK")
 plot!(ca .* 1000, extract(sim, sys.Ca2CaM_N_CaMK), lab="Ca2CaM_N_CaMK")
-plot!(ca .* 1000, extract(sim, sys.Ca4CaM_CaMK), lab="Ca4CaM_CaMK", legend=:topleft, size=(800, 600))
+plot!(ca .* 1000, extract(sim, sys.Ca4CaM_CaMK), lab="Ca4CaM_CaMK", legend=:topleft, size=(600, 600))
 
 #---
 println("Basal activity without ca is ", sol0[sys.CaMKAct][end])
-plot(ca .* 1000, extract(sim, sys.CaMKAct), label=false, title="Active CaMKII", ylims=(0, 0.25); xopts...)
+plot(ca .* 1000, extract(sim, sys.CaMKAct), label=false, title="Active CaMKII", ylims=(0, 0.5); xopts...)
 
 # ## Least-square fitting of steady state CaMKII activity
 @. model_camk(x, p) = p[1] * hil(x, p[2], p[4]) + p[3]
 xdata = ca
 ydata = extract(sim, sys.CaMKAct)
-p0 = [0.2, 1e-3, 0.019, 2.2]
+p0 = [0.2, 1e-3, 0.01, 2.2]
 lb = [0.0, 1e-9, 0.0, 1.0]
-fit = curve_fit(model_camk, xdata, ydata, p0; lower=lb, autodiff=:forwarddiff)
+fit = curve_fit(model_camk, xdata, ydata, inv.(ydata), p0; lower=lb, autodiff=:forwarddiff)
 
 # Parameters
 pestim = coef(fit)
@@ -104,5 +105,5 @@ rmse(fit)
 # Fit result
 yestim = model.(xdata, Ref(pestim))
 p1 = plot(xdata .* 1000, [ydata yestim], lab=["Full model" "Fitted"], line=[:dash :dot], title="CaMKII activity"; xopts...)
-p2 = plot(xdata .* 1000, (yestim .- ydata) ./ ydata .* 100, xlabel="Ca (μM)", title="Relative error (%)", lab=false, yticks=-7:7; xopts...)
+p2 = plot(xdata .* 1000, (yestim .- ydata) ./ ydata .* 100, xlabel="Ca (μM)", title="Relative error (%)", lab=false; xopts...)
 plot(p1, p2, layout=(2, 1), size=(600, 600))

docs/camk-simp.jl

@@ -0,0 +1,57 @@
+# # Test: rapid binding of Ca and calmodulin
+# Assuming rapid equilibrium for CaM0 <--> CaM2 <--> Ca4CaM
+using ModelingToolkit
+using OrdinaryDiffEq
+using SteadyStateDiffEq
+using DiffEqCallbacks
+using Plots
+using LsqFit
+using CaMKIIModel
+using CaMKIIModel: μM, hil
+Plots.default(lw=1.5)
+
+#---
+sys = build_neonatal_ecc_sys(simplify=true, reduce_iso=true)
+sys_simp = build_neonatal_ecc_sys(simplify=true, reduce_iso=true, reduce_camk=true)
+
+println("# of state variable in the original model: ", length(unknowns(sys)))
+println("# of state variable in the simplified model: ", length(unknowns(sys_simp)))
+
+#---
+tend = 500.0
+stimstart = 100.0
+stimend = 300.0
+alg = FBDF()
+@unpack Istim = sys
+callback = build_stim_callbacks(Istim, stimend; period=1, starttime=stimstart)
+prob = ODEProblem(sys, [], tend)
+prob_simp = ODEProblem(sys_simp, [], tend)
+
+#---
+@time sol = solve(prob, alg; callback, abstol=1e-6, reltol=1e-6)
+@time sol_simp = solve(prob_simp, alg; callback, abstol=1e-6, reltol=1e-6)
+
+#---
+plot(sol, idxs=sys.vm*1000, tspan=(295, 300), label= "Full model", title="Action potential")
+plot!(sol_simp, idxs=sys.vm*1000, tspan=(295, 300), label= "Simplified model", line=:dash)
+
+# Something is wrong
+plot(sol, idxs=sys.CaMKAct, label= "Full model", title="Active CaMKII")
+plot!(sol_simp, idxs=sys_simp.CaMKAct, label= "Simplified model", line=:dash)
+
+
+# PCaM too high
+@unpack CaM, KCaM, PCaM, CaMK = sys_simp
+@unpack CaM0, CaM2C, CaM2N, CaM4, CaM0_CaMK, CaM2C_CaMK, CaM2N_CaMK, CaM4_CaMK, CaM0_CaMKP, CaM2C_CaMKP, CaM2N_CaMKP, CaM4_CaMKP = sys_simp
+plot(sol_simp, idxs=[CaM0, CaM2C, CaM2N, CaM4], tspan=(100, 105))
+
+@unpack CaM0, Ca2CaM_C, Ca2CaM_N, Ca4CaM = sys
+plot(sol, idxs=[CaM0, Ca2CaM_C, Ca2CaM_N, Ca4CaM], tspan=(100, 105))
+
+#---
+@parameters Ca = 0 ROS=0
+ksys = CaMKIIModel.get_camkii_fast_ca_binding_sys(Ca; ROS)
+
+for eq in equations(ksys)
+    println(eq)
+end

src/camkii_ros.jl

@@ -8,7 +8,6 @@ ROS activation model: Oxidized Calmodulin Kinase II Regulates Conduction Followi
 function get_camkii_sys(Ca=0μM;
     ROS=0μM,
     binding_To_PCaMK=0,   ## 0.1 for T287D mutation
-    phospho_rate=30Hz,
     name=:camkii_sys,
     simplify=false
 )
@@ -23,7 +22,7 @@ function get_camkii_sys(Ca=0μM;
         k_1N_on = 100Hz / μM    ## 25-260uM-1s-1
         k_1N_off = 2000Hz       ## 1000-4000 s-1
         k_2N_on = 200Hz / μM    ## 50-300uM-1s-1.
-        k_2N_off = 500Hz        ## 500->1000.s-1
+        k_2N_off = 500Hz        ## 500-1000.s-1
 
         ## Ca2+ binding to CaM-CAMKII (KCaM)
         ## C-lobe
@@ -39,13 +38,13 @@ function get_camkii_sys(Ca=0μM;
 
         ## CaM binding to CaMKII
         kCaM0_on = 3.8Hz / mM
-        kCaM2C_on = 0.92Hz / mM
-        kCaM2N_on = 0.12Hz / mM
-        kCaM4_on = 30Hz / mM
         kCaM0_off = 5.5Hz
+        kCaM2C_on = 0.92Hz / μM
         kCaM2C_off = 6.8Hz
+        kCaM2N_on = 0.12Hz / μM
         kCaM2N_off = 1.7Hz
-        kCaM4_off = 1.5Hz
+        kCaM4_on = 30Hz / μM  # 14-60 uM-1s-1
+        kCaM4_off = 1.5Hz  # 1.1 - 2.3 s-1
         kCaM0P_on = kCaM0_on * binding_To_PCaMK
         kCaM2CP_on = kCaM2C_on * binding_To_PCaMK
         kCaM2NP_on = kCaM2N_on * binding_To_PCaMK
@@ -54,12 +53,12 @@ function get_camkii_sys(Ca=0μM;
         kCaM2CP_off = inv(3second)
         kCaM2NP_off = inv(3second)
         kCaM4P_off = inv(3second)
-        k_phosCaM = phospho_rate # 12.6Hz
+        k_phosCaM = 12.6Hz # 30Hz
         k_dephospho = inv(6second)
         k_P1_P2 = inv(60second)
         k_P2_P1 = inv(15second)
 
-        ## Oxidation / reduction
+        ## Oxidation / reduction of Met
         k_BOX = 291Hz / mM
         k_POXP = 291Hz / mM
         k_OXB = inv(45second)
@@ -167,12 +166,10 @@ function get_camkii_sys(Ca=0μM;
     return sys
 end
 
-
 "CaMKII system with ROS activation and fast calcium binding"
 function get_camkii_fast_ca_binding_sys(Ca=0μM;
     ROS=0μM,
     binding_To_PCaMK=0,
-    phospho_rate=30Hz,
     name=:camkii_sys,
     simplify=false
 )
@@ -204,9 +201,9 @@ function get_camkii_fast_ca_binding_sys(Ca=0μM;
 
         ## CaM binding to CaMKII
         kCaM0_on = 3.8Hz / mM
-        kCaM2C_on = 0.92Hz / mM
-        kCaM2N_on = 0.12Hz / mM
-        kCaM4_on = 30Hz / mM
+        kCaM2C_on = 0.92Hz / μM
+        kCaM2N_on = 0.12Hz / μM
+        kCaM4_on = 30Hz / μM
         kCaM0_off = 5.5Hz
         kCaM2C_off = 6.8Hz
         kCaM2N_off = 1.7Hz
@@ -219,7 +216,7 @@ function get_camkii_fast_ca_binding_sys(Ca=0μM;
         kCaM2CP_off = inv(3second)
         kCaM2NP_off = inv(3second)
         kCaM4P_off = inv(3second)
-        k_phosCaM = phospho_rate # 30Hz/12.6Hz
+        k_phosCaM = 12.6Hz # 30Hz
         k_dephospho = inv(6second)
         k_P1_P2 = inv(60second)
         k_P2_P1 = inv(15second)

src/ecc_neonatal.jl

@@ -27,6 +27,7 @@ function build_neonatal_ecc_sys(;
     name=:neonataleccsys,
     simplify=true,
     reduce_iso=false,
+    reduce_camk=false,
 )
     @parameters begin
         ca_o = 1.796mM
@@ -58,7 +59,7 @@ function build_neonatal_ecc_sys(;
     @unpack LCCa_PKAp, LCCb_PKAp, fracPLBp, TnI_PKAp, IKUR_PKAp = barsys
     capdesys = get_ca_pde_sys(; JCa_SR, JCa_SL, TnI_PKAp, rSR_true, rSL_true, dx)
     @unpack Cai_sub_SL, Cai_sub_SR, Cai_mean = capdesys
-    camkiisys = get_camkii_sys(Cai_mean; ROS)
+    camkiisys = reduce_camk ? get_camkii_fast_ca_binding_sys(Cai_mean; ROS) : get_camkii_sys(Cai_mean; ROS)
     icasys = get_ica_sys(na_i, Cai_sub_SL, na_o, ca_o, vm; LCCb_PKAp)
     @unpack INaCa, ICaL, ICaT, ICab = icasys
     inasys = get_ina_sys(na_i, na_o, vm)