use SteadyStateDiffEq

docs/camk-sens.jl

@@ -1,6 +1,7 @@
 # # Sensitivity of CaMKII system to calcium
 using ModelingToolkit
 using OrdinaryDiffEq
+using SteadyStateDiffEq
 using DiffEqCallbacks
 using Plots
 using LsqFit
@@ -15,21 +16,19 @@ sys = get_camkii_sys(Ca; ROS, phospho_rate=0, simplify=true)
 equations(sys)
 
 #---
-tend = 10000.0
-prob = ODEProblem(sys, [], tend)
-alg = Rodas5()
-callback = TerminateSteadyState()
+prob = SteadyStateProblem(sys, [])
+alg = DynamicSS(Rodas5P())
 ca = exp10.(range(log10(0.1μM), log10(100μ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; callback) ## warmup
-sim = solve(EnsembleProblem(prob; prob_func, safetycopy=false), alg; trajectories, callback)
+sol0 = solve(remake(prob, p=[Ca => 0.0]), alg) ## warmup
+sim = solve(EnsembleProblem(prob; prob_func, safetycopy=false), alg; trajectories)
 
 #---
 """Extract values from ensemble simulations by a symbol"""
-extract(sim, k) = map(s -> s[k][end], sim)
+extract(sim, k) = map(s -> s[k], sim)
 """Calculate Root Mean Square Error (RMSE)"""
 rmse(fit) = sqrt(sum(abs2, fit.resid) / length(fit.resid))
 
@@ -47,7 +46,7 @@ plot!(ca .* 1000, extract(sim, sys.Ca4CaM_CaMK), lab="Ca4CaM_CaMK", legend=:topl
 
 #---
 println("Basal activity without ca is ", sol0[sys.CaMKAct][end])
-plot(ca .* 1000, extract(sim, sys.CaMKAct), label=false, title="Active CaMKII"; xopts...)
+plot(ca .* 1000, extract(sim, sys.CaMKAct), label=false, title="Active CaMKII", ylims=(0, 0.25); xopts...)
 
 # ## Least-square fitting of steady state CaMKII activity
 @. model_camk(x, p) = p[1] * hil(x, p[2], p[4]) + p[3]

pluto-nb/test-lcc.jl

@@ -17,7 +17,7 @@ begin
 	using Plots
 	using LsqFit
 	using CaMKIIModel
-	using CaMKIIModel: μM, Hz
+	using CaMKIIModel: μM, Hz, μAμF
 	Plots.default(lw=1.5, fmt=:png)
 end
 
@@ -35,6 +35,7 @@ ps = [
 	sys.ktrCaSR => 50Hz,
 	sys.kRyR => 20Hz,
 	sys.GCaL => 1.3e-4, # 6.3e-5
+	sys.kNaCa => 2.268e-16μAμF / μM^4 # 2.268e-16
 ]
 
 # ╔═╡ 78fb22d4-2f06-49b0-a818-377eedf7ca8d
@@ -47,7 +48,15 @@ prob = ODEProblem(sys, [], tend, ps)
 callback = build_stim_callbacks(Istim, tend; period=1)
 
 # ╔═╡ 761a0baf-b5ec-4caa-ae56-5808a2840987
-sol = solve(prob, FBDF(); callback, abstol=1e-6, reltol=1e-6, progress=true)
+sol = solve(prob, FBDF(); callback, abstol=1e-6, reltol=1e-6, progress=true);
+
+# ╔═╡ 58073dd2-9670-48d1-876e-82ba2a1f5611
+begin
+	# ICaL and INaCa are OK
+	# RyR flux is tiny in NRVM
+	@unpack INaCa, ICaL = sys
+	plot(sol, idxs=[INaCa, ICaL], tspan=(90, 92), ylabel="uA/uF or uM/ms")
+end
 
 # ╔═╡ 8a2eb462-92ee-4ed3-a315-3c58a9a194f4
 plot(sol, idxs=sys.vm*1000, lab="Membrane potential", tspan=(90, 92))
@@ -60,16 +69,8 @@ plot(sol, idxs=[sys.Cai_sub_SL*1E6, sys.Cai_sub_SR*1E6], tspan=(90, 92), ylabel=
 # SR Ca release is reletively small (see JSR)
 plot(sol, idxs=[sys.CaJSR, sys.CaNSR], tspan=(90, 92), ylabel="mM")
 
-# ╔═╡ 58073dd2-9670-48d1-876e-82ba2a1f5611
-begin
-	# ICaL and INaCa are OK
-	# RyR flux is tiny in NRVM
-	@unpack INaCa, ICaL = sys
-	plot(sol, idxs=[INaCa, ICaL], tspan=(90, 92), ylabel="uA/uF or uM/ms")
-end
-
 # ╔═╡ b563ef4c-1dc3-4b3e-9f58-bace24d03d9e
-plot(sol, idxs=sys.CaMKAct, ylim=(0.18, 0.28))
+plot(sol, idxs=sys.CaMKAct, ylim=(0, 0.28))
 
 # ╔═╡ Cell order:
 # ╠═fe7e7c0b-1ae0-4a04-9549-23aba60efcf5
@@ -81,8 +82,8 @@ plot(sol, idxs=sys.CaMKAct, ylim=(0.18, 0.28))
 # ╠═0d97c77a-fdd9-463c-9c53-d40be74b266b
 # ╠═12ded313-721a-4d7b-a2e3-25b313c1a582
 # ╠═761a0baf-b5ec-4caa-ae56-5808a2840987
+# ╠═58073dd2-9670-48d1-876e-82ba2a1f5611
 # ╠═8a2eb462-92ee-4ed3-a315-3c58a9a194f4
 # ╠═152536f6-3cb9-48dd-9467-f1a2b64b3dba
 # ╠═cd11a3d5-bdd8-415c-8b8d-58c4228f8877
-# ╠═58073dd2-9670-48d1-876e-82ba2a1f5611
 # ╠═b563ef4c-1dc3-4b3e-9f58-bace24d03d9e