diff --git a/Manifest.toml b/Manifest.toml index 0f6c53a..da2bfc6 100644 --- a/Manifest.toml +++ b/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" diff --git a/docs/bcl.jl b/docs/bcl.jl index faab5c2..688c3b7 100644 --- a/docs/bcl.jl +++ b/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 diff --git a/docs/camk-sens.jl b/docs/camk-sens.jl index 0a4ecdd..5dbb58f 100644 --- a/docs/camk-sens.jl +++ b/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)) diff --git a/docs/camk-simp.jl b/docs/camk-simp.jl new file mode 100644 index 0000000..3db5e9b --- /dev/null +++ b/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 diff --git a/src/camkii_ros.jl b/src/camkii_ros.jl index 557de2d..85245a1 100644 --- a/src/camkii_ros.jl +++ b/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) diff --git a/src/ecc_neonatal.jl b/src/ecc_neonatal.jl index eb0abe9..af8a7ed 100644 --- a/src/ecc_neonatal.jl +++ b/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)