Coffee (#38)

* pkg * 3Hz * pkg * pkg * Caffeine factor * pkgs * test caffeine * rename * pkg * rename * title * Revert "pkg" This reverts commit 352db14. * wip * pkg * self-hosted * caffeine-effects
ntumitolab · Nov 21, 2024 · 929d2ca · 929d2ca
1 parent a2feff2
commit 929d2ca
Show file tree

Hide file tree

Showing 9 changed files with 144 additions and 117 deletions.
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
@@ -12,8 +12,8 @@ concurrency:
   cancel-in-progress: true
 
 env:
-  NBCONVERT_JOBS: '3'
-  LITERATE_PROC: '3'
+  NBCONVERT_JOBS: '6'
+  LITERATE_PROC: '6'
   JULIA_NUM_THREADS: '2'
   ALLOWERRORS: 'false'
   NBCACHE: '.cache'
@@ -22,7 +22,7 @@ env:
 
 jobs:
   CI:
-    runs-on: ubuntu-latest
+    runs-on: self-hosted
     steps:
     - name: Checkout repository
       uses: actions/checkout@v4

diff --git a/Manifest.toml b/Manifest.toml
diff --git a/docs/_toc.yml b/docs/_toc.yml
@@ -8,6 +8,7 @@ chapters:
   - file: iso-sens
   - file: camk-sens
   - file: camk-simp
-  - file: bcl
-  - file: iso
-  - file: ros
+  - file: effects-bcl
+  - file: iso-effects
+  - file: ros-effects
+  - file: caffeine-effects
diff --git a/docs/caffeine-effects.jl b/docs/caffeine-effects.jl
@@ -0,0 +1,53 @@
+# # Caffeine Effects
+# Caffeine increase RyR opening sensitivity to luminal and subspace calcium
+# In this model, we decrease the mid saturation sub-SR calcium concentration for the opening rate
+using ModelingToolkit
+using OrdinaryDiffEq
+using DiffEqCallbacks
+using Plots
+using CaMKIIModel
+Plots.default(lw=1.5)
+
+#---
+sys = build_neonatal_ecc_sys(simplify=true, reduce_iso=true)
+tend = 500.0
+prob = ODEProblem(sys, [], tend)
+stimstart = 100.0
+stimend = 300.0
+alg = FBDF()
+function add_coffee_affect!(integrator)
+    integrator.ps[sys.RyRsensitivity] = 10
+end
+
+@unpack Istim = sys
+callback = build_stim_callbacks(Istim, stimend; period=1, starttime=stimstart)
+# Add caffeine at t = 200
+callback_caf = CallbackSet(build_stim_callbacks(Istim, stimend; period=1, starttime=stimstart), PresetTimeCallback(200.0, add_coffee_affect!))
+
+#---
+sol = solve(prob, alg; callback, abstol=1e-6, reltol=1e-6)
+sol_caf = solve(prob, alg; callback = callback_caf, abstol=1e-6, reltol=1e-6)
+
+#---
+plot(sol, idxs=sys.vm*1000, title="Action potential", lab="Ctl")
+plot!(sol_caf, idxs=sys.vm*1000, lab="Caf", tspan=(199, 210), ylabel="Voltage (mV)")
+
+#---
+plot(sol, idxs=sys.Cai_sub_SR*1E6, title="Calcium transient (During caffeine addition)", lab="Ctl")
+plot!(sol_caf, idxs=sys.Cai_sub_SR*1E6, tspan=(195, 210), lab="Caf", ylabel="Subspace calcium (nM)")
+
+#---
+plot(sol, idxs=sys.Cai_sub_SR*1E6, title="Calcium transient (After caffeine addition)", lab="Ctl", ylabel="Subspace calcium (nM)")
+plot!(sol_caf, idxs=sys.Cai_sub_SR*1E6, tspan=(295, 300), lab="Caf")
+
+#---
+plot(sol, idxs=sys.CaJSR*1E3, title="SR Calcium (During caffeine addition)", lab="Ctl", ylabel="SR calcium (μM)")
+plot!(sol_caf, idxs=sys.CaJSR*1E3, tspan=(195, 210), lab="Caf", ylims=(200, 850))
+
+#---
+plot(sol, idxs=sys.Jrel, title="RyR Ca flux", lab="Ctl")
+plot!(sol_caf, idxs=sys.Jrel, tspan=(195, 210), lab="Caf", ylabel="μM/ms")
+
+#---
+plot(sol, idxs=sys.CaMKAct, title="Active CaMKII", lab="Ctl")
+plot!(sol_caf, idxs=sys.CaMKAct, lab="Caf")
diff --git a/docs/bcl.jl → docs/effects-bcl.jl b/docs/bcl.jl → docs/effects-bcl.jl
@@ -1,4 +1,4 @@
-# # Pacing frequency
+# # Pacing frequency Effects
 using ModelingToolkit
 using OrdinaryDiffEq
 using Plots
@@ -49,13 +49,32 @@ 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")
 
-# ## Comparing 1-2 Hz
+# ## 3Hz
+# AP not stable
+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=(298, 300))
+
+#---
+plot(sol3, idxs=[sys.Cai_sub_SR, sys.Cai_sub_SL, sys.Cai_mean], tspan=(298, 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", tspan=(299, 300), xlabel="Time (sec.)", ylabel="Voltage (mV)")
+plot!(sol3, idxs=sys.vm*1000, lab="3Hz")
 
 #---
 plot(sol, idxs=sys.CaMKAct, title="CaMKII", lab="1Hz")
 plot!(sol2, idxs=sys.CaMKAct, lab="2Hz",  ylim=(0.0, 1.0), xlabel="Time (sec.)", ylabel="Act. fraction (AU)")
+plot!(sol3, idxs=sys.CaMKAct, lab="3Hz")
 
 # ## Data fitting
 ### Pacing duration and CaMKII activity

diff --git a/docs/iso.jl → docs/iso-effects.jl b/docs/iso.jl → docs/iso-effects.jl
diff --git a/docs/ros.jl → docs/ros-effects.jl b/docs/ros.jl → docs/ros-effects.jl
diff --git a/pluto-nb/test-lcc.jl b/pluto-nb/test-lcc.jl
@@ -15,7 +15,6 @@ begin
 	using OrdinaryDiffEq
 	using DiffEqCallbacks
 	using Plots
-	using LsqFit
 	using CaMKIIModel
 	using CaMKIIModel: μM, Hz, μAμF
 	Plots.default(lw=1.5)
@@ -28,14 +27,11 @@ Base.current_project()
 sys = build_neonatal_ecc_sys(simplify=true, reduce_iso=true);
 
 # ╔═╡ fb697567-fa32-4963-a490-0668657b22f5
-tend = 100.0
+tend = 500.0
 
 # ╔═╡ 638be248-21ad-4833-bbc8-f2f9a805a900
 ps = [
-	sys.ktrCaSR => 50Hz,
-	sys.kRyR => 20Hz,
-	sys.GCaL => 1e-4, # 6.3e-5
-	sys.kNaCa => 2.268e-16μAμF / μM^4 # 2.268e-16
+	sys.RyRsensitivity => 10
 ]
 
 # ╔═╡ 78fb22d4-2f06-49b0-a818-377eedf7ca8d
@@ -58,6 +54,9 @@ begin
 	plot(sol, idxs=[INaCa, ICaL], tspan=(90, 92), ylabel="uA/uF or uM/ms")
 end
 
+# ╔═╡ 32ea73cb-f2f1-4530-a401-f7db9fb1b313
+plot(sol, idxs=sys.PO1RyR, tspan=(90, 92))
+
 # ╔═╡ 8a2eb462-92ee-4ed3-a315-3c58a9a194f4
 plot(sol, idxs=sys.vm*1000, lab="Membrane potential", tspan=(90, 92))
 
@@ -66,28 +65,11 @@ plot(sol, idxs=sys.vm*1000, lab="Membrane potential", tspan=(90, 92))
 plot(sol, idxs=[sys.Cai_sub_SL*1E6, sys.Cai_sub_SR*1E6], tspan=(90, 92), ylabel="nM", lab=["Ca (sub SL)" "Ca (sub SR)"], ylims=(100, 800))
 
 # ╔═╡ cd11a3d5-bdd8-415c-8b8d-58c4228f8877
-# SR Ca release is reletively small (see JSR)
 plot(sol, idxs=[sys.CaJSR, sys.CaNSR], tspan=(90, 92), ylabel="mM")
 
 # ╔═╡ b563ef4c-1dc3-4b3e-9f58-bace24d03d9e
 plot(sol, idxs=sys.CaMKAct)
 
-# ╔═╡ 3e63b27f-a13a-45a6-9f95-83d31d64db34
-begin 
-	@unpack CaM0_CaMK, Ca2CaM_C_CaMK, Ca2CaM_N_CaMK, Ca4CaM_CaMK, CaM0_CaMKP, Ca2CaM_C_CaMKP, Ca2CaM_N_CaMKP, Ca4CaM_CaMKP, CaMKP, CaMKP2, CaMK= sys;
-	KCaM = CaM0_CaMK + Ca2CaM_C_CaMK + Ca2CaM_N_CaMK + Ca4CaM_CaMK
-	PCaM = CaM0_CaMKP + Ca2CaM_C_CaMKP + Ca2CaM_N_CaMKP + Ca4CaM_CaMKP
-end;
-
-# ╔═╡ 5d7497ec-d7b0-41b6-ba68-4015bb98090a
-begin
-	plot(sol, idxs=KCaM, lab="KCaM")
-	plot!(sol, idxs=PCaM, lab="PCaM")
-	plot!(sol, idxs=CaMKP)
-	plot!(sol, idxs=CaMKP2)
-	plot!(sol, idxs=CaMK, lab="Inactive CaMK")
-end
-
 # ╔═╡ Cell order:
 # ╠═fe7e7c0b-1ae0-4a04-9549-23aba60efcf5
 # ╠═f4d83570-9d1c-11ef-0654-53e2d6f12d6e
@@ -99,9 +81,8 @@ end
 # ╠═12ded313-721a-4d7b-a2e3-25b313c1a582
 # ╠═761a0baf-b5ec-4caa-ae56-5808a2840987
 # ╠═58073dd2-9670-48d1-876e-82ba2a1f5611
+# ╠═32ea73cb-f2f1-4530-a401-f7db9fb1b313
 # ╠═8a2eb462-92ee-4ed3-a315-3c58a9a194f4
 # ╠═152536f6-3cb9-48dd-9467-f1a2b64b3dba
 # ╠═cd11a3d5-bdd8-415c-8b8d-58c4228f8877
 # ╠═b563ef4c-1dc3-4b3e-9f58-bace24d03d9e
-# ╠═3e63b27f-a13a-45a6-9f95-83d31d64db34
-# ╠═5d7497ec-d7b0-41b6-ba68-4015bb98090a
diff --git a/src/er.jl b/src/er.jl
@@ -5,16 +5,17 @@ function get_ser_sys(Cai_sub_SR; fracPLB_CKp=0, fracPLBp=0, RyR_CKp=0.2, V_sub_S
         VNSR = 0.9 * VSR
         VJSR = VSR - VNSR
         # RyR
-        kRyR = 20Hz
+        kRyR = 20Hz # 10 Hz
         kaposRyR = 1000Hz
         kanegRyR = 160Hz
+        RyRsensitivity = 1.0
         # SERCA
         VmaxSR = 0.9996mM*Hz
         KmfSR = 0.5μM
         KmrSR = 3.5mM
         kSRleak = 0.005Hz
         fracPKA_PLBo = 1 - 0.079755
-        ktrCaSR = 50Hz
+        ktrCaSR = 50Hz  # 5Hz
         csqntot = 24.750mM
         Kmcsqn = 0.8mM
     end
@@ -36,7 +37,7 @@ function get_ser_sys(Cai_sub_SR; fracPLB_CKp=0, fracPLBp=0, RyR_CKp=0.2, V_sub_S
     fCKII_PLB = (1 - 0.5 * fracPLB_CKp)  # Max effect: fCKII_PLB=0.5
     fPKA_PLB = ((1 - fracPLBp) / fracPKA_PLBo) * (1 - 0.5531) + 0.5531
     # Select the smaller value (resulting in max reduction of Kmf)
-    Kmfp = KmfSR * min(fCKII_PLB, fPKA_PLB)  #fCKII_PLB
+    Kmfp = KmfSR * min(fCKII_PLB, fPKA_PLB)  # fCKII_PLB
     fSR = (Cai_sub_SR / Kmfp)^2
     rSR = (CaNSR / KmrSR)^2
     kleak = (1 + 5 * RyR_CKp) * kSRleak / 2
@@ -45,7 +46,7 @@ function get_ser_sys(Cai_sub_SR; fracPLB_CKp=0, fracPLBp=0, RyR_CKp=0.2, V_sub_S
         1 ~ PO1RyR + PC1RyR,
         Jrel ~ kRyR * PO1RyR * (CaJSR - Cai_sub_SR),
         KmRyR ~ (3.51 / (1 + exp((CaJSR - 530μM) / 200μM)) + 0.25) * μM,
-        D(PO1RyR) ~ kaposRyR * hil(Cai_sub_SR, KmRyR, 4) * PC1RyR - kanegRyR * PO1RyR,
+        D(PO1RyR) ~ kaposRyR * hil(Cai_sub_SR, KmRyR / RyRsensitivity, 4) * PC1RyR - kanegRyR * PO1RyR,
         Jup ~ VmaxSR * (fSR - rSR) / (1 + fSR + rSR),
         Jleak ~ kleak * (CaNSR - Cai_sub_SR),
         Jtr ~ ktrCaSR * (CaNSR - CaJSR),