diff --git a/Manifest.toml b/Manifest.toml
index 1b14fc1..d566acf 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -5,9 +5,9 @@ manifest_format = "2.0"
 project_hash = "f9d4542ded0a8d55955a6ceb5bc260f4265a0a2c"
 
 [[deps.ADTypes]]
-git-tree-sha1 = "5a5eafb8344b81b8c2237f8a6f6b3602b3f6180e"
+git-tree-sha1 = "eea5d80188827b35333801ef97a40c2ed653b081"
 uuid = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
-version = "1.8.1"
+version = "1.9.0"
 weakdeps = ["ChainRulesCore", "EnzymeCore"]
 
     [deps.ADTypes.extensions]
@@ -510,9 +510,9 @@ uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 
 [[deps.Distributions]]
 deps = ["AliasTables", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsAPI", "StatsBase", "StatsFuns"]
-git-tree-sha1 = "e6c693a0e4394f8fda0e51a5bdf5aef26f8235e9"
+git-tree-sha1 = "d7477ecdafb813ddee2ae727afa94e9dcb5f3fb0"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
-version = "0.25.111"
+version = "0.25.112"
 
     [deps.Distributions.extensions]
     DistributionsChainRulesCoreExt = "ChainRulesCore"
@@ -555,9 +555,9 @@ version = "0.6.0"
 
 [[deps.DynamicQuantities]]
 deps = ["Compat", "DispatchDoctor", "PackageExtensionCompat", "TestItems", "Tricks"]
-git-tree-sha1 = "317430fb3dcfb17831eb0e4bcd852358a48a4e7a"
+git-tree-sha1 = "82b3ab10e1d14cd79d813d30005d52ac3ef9ea97"
 uuid = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
-version = "1.0.0"
+version = "1.1.0"
 
     [deps.DynamicQuantities.extensions]
     DynamicQuantitiesLinearAlgebraExt = "LinearAlgebra"
@@ -622,9 +622,9 @@ version = "0.8.5"
 
 [[deps.FFMPEG]]
 deps = ["FFMPEG_jll"]
-git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8"
+git-tree-sha1 = "53ebe7511fa11d33bec688a9178fac4e49eeee00"
 uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
-version = "0.4.1"
+version = "0.4.2"
 
 [[deps.FFMPEG_jll]]
 deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "PCRE2_jll", "Zlib_jll", "libaom_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
@@ -915,9 +915,9 @@ version = "0.21.4"
 
 [[deps.JpegTurbo_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "c84a835e1a09b289ffcd2271bf2a337bbdda6637"
+git-tree-sha1 = "25ee0be4d43d0269027024d75a24c24d6c6e590c"
 uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
-version = "3.0.3+0"
+version = "3.0.4+0"
 
 [[deps.JuliaFormatter]]
 deps = ["CSTParser", "CommonMark", "DataStructures", "Glob", "PrecompileTools", "TOML", "Tokenize"]
@@ -964,9 +964,9 @@ version = "18.1.7+0"
 
 [[deps.LZO_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "70c5da094887fd2cae843b8db33920bac4b6f07d"
+git-tree-sha1 = "854a9c268c43b77b0a27f22d7fab8d33cdb3a731"
 uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
-version = "2.10.2+0"
+version = "2.10.2+1"
 
 [[deps.LaTeXStrings]]
 git-tree-sha1 = "50901ebc375ed41dbf8058da26f9de442febbbec"
@@ -1265,9 +1265,9 @@ uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 
 [[deps.ModelingToolkit]]
 deps = ["AbstractTrees", "ArrayInterface", "BlockArrays", "Combinatorics", "Compat", "ConstructionBase", "DataStructures", "DiffEqBase", "DiffEqCallbacks", "DiffEqNoiseProcess", "DiffRules", "Distributed", "Distributions", "DocStringExtensions", "DomainSets", "DynamicQuantities", "ExprTools", "Expronicon", "FindFirstFunctions", "ForwardDiff", "FunctionWrappers", "FunctionWrappersWrappers", "Graphs", "InteractiveUtils", "JuliaFormatter", "JumpProcesses", "Latexify", "Libdl", "LinearAlgebra", "MLStyle", "NaNMath", "NonlinearSolve", "OrderedCollections", "PrecompileTools", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLBase", "SciMLStructures", "Serialization", "Setfield", "SimpleNonlinearSolve", "SparseArrays", "SpecialFunctions", "StaticArrays", "SymbolicIndexingInterface", "SymbolicUtils", "Symbolics", "URIs", "UnPack", "Unitful"]
-git-tree-sha1 = "ebcdaf39b0cee8051193911fe2c97f327fa6735f"
+git-tree-sha1 = "00140c3b93fbbbd317a71a65c578306d0d72b1be"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
-version = "9.40.1"
+version = "9.41.0"
 
     [deps.ModelingToolkit.extensions]
     MTKBifurcationKitExt = "BifurcationKit"
@@ -1922,9 +1922,9 @@ version = "0.6.43"
 
 [[deps.SciMLBase]]
 deps = ["ADTypes", "Accessors", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "Expronicon", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "PrecompileTools", "Preferences", "Printf", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLOperators", "SciMLStructures", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables"]
-git-tree-sha1 = "82584ea03bda16156ceaa6af75f9c8e1287029d3"
+git-tree-sha1 = "71857d6bab17e7ac6802d86ffcc75423b8c1d812"
 uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
-version = "2.53.2"
+version = "2.54.0"
 
     [deps.SciMLBase.extensions]
     SciMLBaseChainRulesCoreExt = "ChainRulesCore"
@@ -1994,9 +1994,9 @@ version = "1.2.0"
 
 [[deps.SimpleNonlinearSolve]]
 deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "DiffEqBase", "DiffResults", "DifferentiationInterface", "FastClosures", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "MaybeInplace", "PrecompileTools", "Reexport", "SciMLBase", "Setfield", "StaticArraysCore"]
-git-tree-sha1 = "536c0ee0b4b766ddee24220c6bb60932df4e2c39"
+git-tree-sha1 = "97a7e10717bd27e600ba55da9d90a5628c1a403b"
 uuid = "727e6d20-b764-4bd8-a329-72de5adea6c7"
-version = "1.12.1"
+version = "1.12.2"
 
     [deps.SimpleNonlinearSolve.extensions]
     SimpleNonlinearSolveChainRulesCoreExt = "ChainRulesCore"
@@ -2057,9 +2057,9 @@ version = "2.20.0"
 
 [[deps.SparseMatrixColorings]]
 deps = ["ADTypes", "Compat", "DataStructures", "DocStringExtensions", "LinearAlgebra", "Random", "SparseArrays"]
-git-tree-sha1 = "996dff77d814c45c3f2342fa0113e4ad31e712e8"
+git-tree-sha1 = "1b17167b4a50197a6e17192acb233e82941b791f"
 uuid = "0a514795-09f3-496d-8182-132a7b665d35"
-version = "0.4.0"
+version = "0.4.1"
 
 [[deps.Sparspak]]
 deps = ["Libdl", "LinearAlgebra", "Logging", "OffsetArrays", "Printf", "SparseArrays", "Test"]
@@ -2179,9 +2179,9 @@ version = "5.2.2+0"
 
 [[deps.SymbolicIndexingInterface]]
 deps = ["Accessors", "ArrayInterface", "RuntimeGeneratedFunctions", "StaticArraysCore"]
-git-tree-sha1 = "988e04b34a4c3b824fb656f542473df99a4f610d"
+git-tree-sha1 = "0225f7c62f5f78db35aae6abb2e5cabe38ce578f"
 uuid = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
-version = "0.3.30"
+version = "0.3.31"
 
 [[deps.SymbolicLimits]]
 deps = ["SymbolicUtils"]
@@ -2205,9 +2205,9 @@ version = "3.7.1"
 
 [[deps.Symbolics]]
 deps = ["ADTypes", "ArrayInterface", "Bijections", "CommonWorldInvalidations", "ConstructionBase", "DataStructures", "DiffRules", "Distributions", "DocStringExtensions", "DomainSets", "DynamicPolynomials", "IfElse", "LaTeXStrings", "LambertW", "Latexify", "Libdl", "LinearAlgebra", "LogExpFunctions", "MacroTools", "Markdown", "NaNMath", "PrecompileTools", "Primes", "RecipesBase", "Reexport", "RuntimeGeneratedFunctions", "SciMLBase", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArraysCore", "SymbolicIndexingInterface", "SymbolicLimits", "SymbolicUtils", "TermInterface"]
-git-tree-sha1 = "8b48697e7fec6d4b7c4a9fe892857a5ed2bae7e8"
+git-tree-sha1 = "e0af193e66eb5c1767f717437c0131d51f539f8a"
 uuid = "0c5d862f-8b57-4792-8d23-62f2024744c7"
-version = "6.12.0"
+version = "6.13.0"
 
     [deps.Symbolics.extensions]
     SymbolicsForwardDiffExt = "ForwardDiff"
@@ -2556,9 +2556,9 @@ version = "1.2.13+1"
 
 [[deps.Zstd_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "e678132f07ddb5bfa46857f0d7620fb9be675d3b"
+git-tree-sha1 = "555d1076590a6cc2fdee2ef1469451f872d8b41b"
 uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
-version = "1.5.6+0"
+version = "1.5.6+1"
 
 [[deps.eudev_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "gperf_jll"]
diff --git a/docs/_toc.yml b/docs/_toc.yml
index a70caec..57b0ea9 100644
--- a/docs/_toc.yml
+++ b/docs/_toc.yml
@@ -4,7 +4,8 @@
 format: jb-book
 root: index
 chapters:
-  - file: iso-sense
+  - file: iso-sens
+  - file: camk-sens
   - file: camkii_calwave
   - file: camkii_ros
   - file: bcl
diff --git a/docs/bcl.ipynb b/docs/bcl.ipynb
deleted file mode 100644
index ceafe85..0000000
--- a/docs/bcl.ipynb
+++ /dev/null
@@ -1,226 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Pacing frequency"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 3,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "using ModelingToolkit\n",
-    "using DifferentialEquations\n",
-    "using Plots\n",
-    "using CaMKIIModel\n",
-    "Plots.default(fmt=:png)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "sys = build_neonatal_ecc_sys(simplify=true)\n",
-    "tend = 300.0\n",
-    "prob = ODEProblem(sys, [], tend)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## 1Hz"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "@unpack Istim = sys\n",
-    "callback = build_stim_callbacks(Istim, tend; period=1)\n",
-    "alg = FBDF()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "sol = solve(prob, alg; callback, abstol=1e-6, reltol=1e-6, maxiters=Int(1e8))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol, idxs=sys.vm, tspan=(295, 300), title=\"Action potential\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol, idxs=[sys.Cai_sub_SR, sys.Cai_sub_SL, sys.Cai_mean], tspan=(299, 300), title=\"Calcium transcient\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol, idxs=sys.CaMKAct, title=\"Active CaMKII\")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## 2Hz"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "callback = build_stim_callbacks(Istim, tend; period=1/2)\n",
-    "sol2 = solve(prob, alg; callback, abstol=1e-6, reltol=1e-6, maxiters=Int(1e8))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol2, idxs=sys.vm, tspan=(298, 300), title=\"Action potential\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol2, idxs=[sys.Cai_sub_SR, sys.Cai_sub_SL, sys.Cai_mean], tspan=(299, 300), title=\"Calcium transcient\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol2, idxs=sys.CaMKAct, title=\"Active CaMKII\")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## 3Hz"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "callback = build_stim_callbacks(Istim, tend; period=1/3)\n",
-    "@time sol3 = solve(prob, alg; callback, abstol=1e-6, reltol=1e-6, maxiters=Int(1e8))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol3, idxs=sys.vm, tspan=(298, 300), title=\"Action potential\")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "The Na channel is not fully recovered."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol3, idxs=[sys.i_Nam, sys.i_Nah, sys.i_Naj], tspan=(299, 300), title=\"Sodium gating\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol3, idxs=[sys.Cai_sub_SR, sys.Cai_sub_SL, sys.Cai_mean], tspan=(299, 300), title=\"Calcium transcient\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol3, idxs=sys.CaMKAct, title=\"Active CaMKII\")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Comparing 1-3 Hz"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol, idxs=sys.CaMKAct, title=\"Active CaMKII\", lab=\"1Hz\")\n",
-    "plot!(sol2, idxs=sys.CaMKAct, lab=\"2Hz\")\n",
-    "plot!(sol3, idxs=sys.CaMKAct, lab=\"3Hz\", ylim=(0.0, 0.8))"
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Julia 1.10.5",
-   "language": "julia",
-   "name": "julia-1.10"
-  },
-  "language_info": {
-   "file_extension": ".jl",
-   "mimetype": "application/julia",
-   "name": "julia",
-   "version": "1.10.5"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/docs/bcl.jl b/docs/bcl.jl
new file mode 100644
index 0000000..091aa47
--- /dev/null
+++ b/docs/bcl.jl
@@ -0,0 +1,60 @@
+# # Pacing frequency
+
+using ModelingToolkit
+using DifferentialEquations
+using Plots
+using CaMKIIModel
+
+# Setup system
+sys = build_neonatal_ecc_sys(simplify=true)
+tend = 300.0
+prob = ODEProblem(sys, [], tend)
+
+# ## 1Hz
+@unpack Istim = sys
+callback = build_stim_callbacks(Istim, tend; period=1)
+alg = FBDF()
+sol = solve(prob, alg; callback, abstol=1e-6, reltol=1e-6, maxiters=Int(1e8))
+
+#---
+plot(sol, idxs=sys.vm, tspan=(295, 300), title="Action potential")
+
+#---
+plot(sol, idxs=[sys.Cai_sub_SR, sys.Cai_sub_SL, sys.Cai_mean], tspan=(299, 300), title="Calcium transcient")
+
+#---
+plot(sol, idxs=sys.CaMKAct, title="Active CaMKII")
+
+# ## 2Hz
+callback = build_stim_callbacks(Istim, tend; period=1/2)
+sol2 = solve(prob, alg; callback, abstol=1e-6, reltol=1e-6, maxiters=Int(1e8))
+
+#---
+plot(sol2, idxs=sys.vm, tspan=(298, 300), title="Action potential")
+
+#---
+plot(sol2, idxs=[sys.Cai_sub_SR, sys.Cai_sub_SL, sys.Cai_mean], tspan=(299, 300), title="Calcium transcient")
+
+#---
+plot(sol2, idxs=sys.CaMKAct, title="Active CaMKII")
+
+# ## 3Hz
+callback = build_stim_callbacks(Istim, tend; period=1/3)
+@time sol3 = solve(prob, alg; callback, abstol=1e-6, reltol=1e-6, maxiters=Int(1e8))
+
+#---
+plot(sol3, idxs=sys.vm, tspan=(298, 300), title="Action potential")
+
+# The Na channel is not fully recovered.
+plot(sol3, idxs=[sys.i_Nam, sys.i_Nah, sys.i_Naj], tspan=(299, 300), title="Sodium gating")
+
+#---
+plot(sol3, idxs=[sys.Cai_sub_SR, sys.Cai_sub_SL, sys.Cai_mean], tspan=(299, 300), title="Calcium transcient")
+
+#---
+plot(sol3, idxs=sys.CaMKAct, title="Active CaMKII")
+
+# ## Comparing 1-3 Hz
+plot(sol, idxs=sys.CaMKAct, title="Active CaMKII", lab="1Hz")
+plot!(sol2, idxs=sys.CaMKAct, lab="2Hz")
+plot!(sol3, idxs=sys.CaMKAct, lab="3Hz", ylim=(0.0, 0.8))
diff --git a/docs/camk-sens.jl b/docs/camk-sens.jl
new file mode 100644
index 0000000..2ca2a4b
--- /dev/null
+++ b/docs/camk-sens.jl
@@ -0,0 +1,63 @@
+# # Sensitivity of CaMKII system to calcium
+using ModelingToolkit
+using DifferentialEquations
+using Plots
+using LsqFit
+using CaMKIIModel
+using CaMKIIModel: μM, hil
+
+# ## 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)
+
+#---
+tend=100.0
+prob = ODEProblem(sys, [], tend)
+alg = Rodas5()
+callback = TerminateSteadyState()
+ca = exp10.(range(log10(1e-2μM), log10(10μM), length=1001))
+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)
+
+"""Extract values from ensemble simulations by a symbol"""
+extract(sim, k) = map(s->s[k][end], sim)
+
+
+# Components
+plot(ca .* 1000, extract(sim, sys.CaM0), lab="CaM0", ylabel="Conc. (mM)", xlabel="Ca (μM)", xscale=:log10)
+plot!(ca .* 1000, extract(sim, sys.Ca2CaM_C), lab="Ca2CaM_C")
+plot!(ca .* 1000, extract(sim, sys.Ca2CaM_N), lab="Ca2CaM_N")
+plot!(ca .* 1000, extract(sim, sys.Ca4CaM), lab="Ca4CaM")
+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, 800))
+
+# Active CaMKII
+println("Basal activity when Ca=0 is ", sol0[sys.CaMKAct][end])
+plot(ca .* 1000, extract(sim, sys.CaMKAct), lab="Active CaMKII", xlabel="Ca (μM)", xscale=:log10)
+
+# ## Least-square fitting of CaMKII activity
+hil_s(x, k, n) = sign(x*k) * (abs(x)^n) / (abs(x)^n + abs(k)^n)
+@. model_camk(x, p) = p[1] * hil_s(x, p[2], p[3]) + p[4]
+xdata = ca
+ydata = extract(sim, sys.CaMKAct)
+
+p0 = [0.4, 1e-3, 3.0, 0.01]
+lb = [0.0, 0.0, 0.1, 0.0]
+
+fit = curve_fit(model_camk, xdata, ydata, p0; lower=lb)
+pestim = coef(fit)
+
+#---
+confidence_inter = confint(fit; level=0.95)
+
+#---
+plot(xdata, ydata, lab="Data", line=:dash, title="CaMKII activity", xscale=:log10)
+plot!(x-> model_camk(x, pestim), xdata, lab="Fitted", line=:dot, legend=:topleft)
diff --git a/docs/iso-sens.jl b/docs/iso-sens.jl
new file mode 100644
index 0000000..fa0c108
--- /dev/null
+++ b/docs/iso-sens.jl
@@ -0,0 +1,118 @@
+# # Sensitivity to ISO
+using ModelingToolkit
+using DifferentialEquations
+using Plots
+using LsqFit
+using CaMKIIModel
+using CaMKIIModel: μM
+Plots.default(lw=2)
+
+# ## Setup b1AR system
+@parameters ATP=5000μM ISO=0μM
+sys = get_bar_sys(ATP, ISO; simplify=true)
+
+#---
+tend=10000.0
+prob = ODEProblem(sys, [], tend)
+alg = Rodas5()
+callback = TerminateSteadyState()
+# Log scale for ISO concentration
+iso = exp10.(range(log10(1e-10μM), log10(1μM), length=1001))
+
+prob_func = (prob, i, repeat) -> begin
+    remake(prob, p=[ISO => iso[i]])
+end
+trajectories = length(iso)
+sol = solve(prob, alg; callback) ## warmup
+sim = solve(EnsembleProblem(prob; prob_func, safetycopy=false), alg; trajectories, callback)
+
+# ## time series
+# It took over 1000 seconds to reach equilibrium.
+plot(sol, idxs=sys.PKACI/sys.RItot)
+
+#---
+"""Extract values from ensemble simulations by a symbol"""
+extract(sim, k) = map(s->s[k][end], sim)
+
+#---
+plot(iso .* 1000, extract(sim, sys.cAMP), lab="cAMP", ylabel="Conc. (mM)", xlabel="ISO (μM)")
+
+#---
+plot(iso .* 1000, extract(sim, sys.PKACI/sys.RItot), lab="PKACI", xlabel="ISO (μM)", ylabel="Activation fraction")
+plot!(iso .* 1000, extract(sim, sys.PKACII/sys.RIItot), lab="PKACII")
+
+# ## Least-square fitting of PKACI activity
+@. model(x, p) = p[1] * x / (x + p[2]) + p[3]
+xdata = iso
+ydata = extract(sim, sys.PKACI/sys.RItot)
+p0 = [0.3, 1E-5, 0.08]
+lb = [0.0, 0.0, 0.0]
+pkac1_fit = curve_fit(model, xdata, ydata, p0; lower=lb)
+pkac1_coef = coef(pkac1_fit)
+
+println("PKACI")
+println("Basal activity: ", pkac1_coef[3])
+println("Maximal activity: ", pkac1_coef[1] + pkac1_coef[3])
+println("Michaelis constant: ", pkac1_coef[2] * 1000, " μM")
+confidence_inter = confint(pkac1_fit; level=0.95)
+
+#---
+plot(xdata, ydata, lab="Data", line=:dash, title="PKACI", xscale=:log10)
+plot!(x-> model(x, pkac1_coef), xdata, lab="Fitted", line=:dot, legend=:topleft)
+
+# ## Least-square fitting of PKACII activity
+xdata = iso
+ydata = extract(sim, sys.PKACII/sys.RIItot)
+p0 = [0.4, 0.01μM, 0.2]
+lb = [0.0, 0.0, 0.0]
+pkac2_fit = curve_fit(model, xdata, ydata, p0; lower=lb)
+pkac2_coef = coef(pkac2_fit)
+
+println("PKACII")
+println("Basal activity: ", pkac2_coef[3])
+println("Maximal activity: ", pkac2_coef[1] + pkac2_coef[3])
+println("Michaelis constant: ", pkac2_coef[2] * 1000, " μM")
+confidence_inter = confint(pkac2_fit; level=0.95)
+
+#---
+plot(xdata, ydata, lab="Data", line=:dash, title="PKACII", xscale=:log10)
+plot!(x-> model(x, pkac2_coef), xdata, lab="Fitted", line=:dot, legend=:topleft)
+
+# ## Least-square fitting of PP1 activity
+@. model_pp1(x, p) = p[1] * p[2] / (x + p[2]) + p[3]
+xdata = iso
+ydata = extract(sim, sys.PP1/sys.PP1tot)
+p0 = [0.1, 8e-6, 0.8]
+lb = [0.0, 0.0, 0.0]
+pp1_fit = curve_fit(model_pp1, xdata, ydata, p0; lower=lb)
+pp1_coef = coef(pp1_fit)
+
+println("PP1")
+println("Basal activity: ", pp1_coef[1] + pp1_coef[3])
+println("Minimal activity: ",  pp1_coef[3])
+println("Repressive Michaelis constant: ", pp1_coef[2] * 1000, " μM")
+confidence_inter = confint(pp1_fit; level=0.95)
+
+#---
+plot(xdata, ydata, lab="Data", line=:dash, title="PP1", xscale=:log10)
+plot!(x-> model_pp1(x, pp1_coef), xdata, lab="Fitted", line=:dot)
+
+# ## Least-square fitting of PLBp
+# Signed Hill function
+hil_s(x, k, n) = sign(x*k) * (abs(x)^n) / (abs(x)^n + abs(k)^n)
+xdata = iso
+ydata = extract(sim, sys.PLBp/sys.PLBtot)
+@. model_plb(x, p) = p[1] * hil_s(x, p[2], p[3]) + p[4]
+p0 = [0.8, 2e-6, 1.0, 0.1]
+lb = [0.0, 0.0, 0.1, -Inf]
+fit = curve_fit(model_plb, xdata, ydata, p0; lower=lb)
+pestim = coef(fit)
+
+#---
+confidence_inter = confint(fit; level=0.95)
+
+# It does not fit as good as the previous ones because it's in fact a cubic equation.
+plot(xdata, ydata, lab="Data", line=:dash, title="PLBp", xscale=:log10)
+plot!(x-> model_plb(x, pestim), xdata, lab="Fitted", line=:dot, legend=:topleft)
+#---
+plot(xdata, (model_plb.(xdata, Ref(pestim)) .- ydata)./ydata, lab="Error", xscale=:log10)
diff --git a/docs/iso-sense.jl b/docs/iso-sense.jl
deleted file mode 100644
index 5325e07..0000000
--- a/docs/iso-sense.jl
+++ /dev/null
@@ -1,113 +0,0 @@
-# # Sensitivity to ISO
-using ModelingToolkit
-using DifferentialEquations
-using Plots
-using LsqFit
-using CaMKIIModel
-using CaMKIIModel: μM, hil
-Plots.default(lw=2)
-
-# ##    Setup b1AR system
-@parameters ATP=5000μM ISO=0μM
-sys = get_bar_sys(ATP, ISO; simplify=true)
-
-#---
-tend=100000.0
-prob = ODEProblem(sys, [], tend)
-alg = Rodas5()
-callback = TerminateSteadyState()
-iso = range(0.0, 0.5μM, length=1001)
-prob_func = (prob, i, repeat) -> begin
-    remake(prob, p=[ISO => iso[i]])
-end
-trajectories = length(iso)
-sol = solve(prob, alg; callback) ## warmup
-sim = solve(EnsembleProblem(prob; prob_func, safetycopy=false), alg; trajectories, callback)
-
-# ## time series
-# It took over 1000 seconds to reach equilibrium.
-plot(sim[1], idxs=sys.PKACI/sys.RItot)
-
-#---
-"""Extract values from ensemble simulations by a symbol"""
-extract(sim, k) = map(s->s[k][end], sim)
-
-#---
-plot(iso .* 1000, extract(sim, sys.cAMP), lab="cAMP", ylabel="Conc. (mM)", xlabel="ISO (μM)")
-
-#---
-plot(iso .* 1000, extract(sim, sys.PKACI/sys.RItot), lab="PKACI", xlabel="ISO (μM)", ylabel="Activation fraction")
-plot!(iso .* 1000, extract(sim, sys.PKACII/sys.RIItot), lab="PKACII")
-
-# ## Least-square fitting of PKACI activity
-@. model(x, p) = p[1] * x / (x + p[2]) + p[3]
-xdata = iso
-ydata = extract(sim, sys.PKACI/sys.RItot)
-p0 = [0.3, 0.04, 0.08]
-lb = [0.0, 0.0, 0.0]
-fit = curve_fit(model, xdata, ydata, p0; lower=lb)
-pestim = coef(fit)
-
-println("PKACI")
-println("Basal activity: ", pestim[3])
-println("Maximal activity: ", pestim[1] + pestim[3])
-println("Michaelis constant: ", pestim[2] * 1000, " μM")
-confidence_inter = confint(fit; level=0.95)
-
-#---
-plot(xdata, ydata, lab="Data", line=:dash, title="PKACI")
-plot!(x-> model(x, pestim), xdata, lab="Fitted", line=:dot)
-
-# ## Least-square fitting of PKACII activity
-xdata = iso
-ydata = extract(sim, sys.PKACII/sys.RIItot)
-p0 = [0.4, 0.04, 0.2]
-fit = curve_fit(model, xdata, ydata, p0; lower=lb)
-pestim = coef(fit)
-
-println("PKACII")
-println("Basal activity: ", pestim[3])
-println("Maximal activity: ", pestim[1] + pestim[3])
-println("Michaelis constant: ", pestim[2] * 1000, " μM")
-confidence_inter = confint(fit; level=0.95)
-
-#---
-plot(xdata, ydata, lab="Data", line=:dash, title="PKACII")
-plot!(x-> model(x, pestim), xdata, lab="Fitted", line=:dot)
-
-# ## Least-square fitting of PP1 activity
-@. model_pp1(x, p) = p[1] * p[2] / (x + p[2]) + p[3]
-xdata = iso
-ydata = extract(sim, sys.PP1/sys.PP1tot)
-p0 = [0.1, 0.00003, 0.8]
-lb = [0.0, 0.0, 0.0]
-fit = curve_fit(model_pp1, xdata, ydata, p0; lower=lb)
-pestim = coef(fit)
-
-println("PP1")
-println("Basal activity: ", pestim[1] + pestim[3])
-println("Minimal activity: ",  pestim[3])
-println("Repressive Michaelis constant: ", pestim[2] * 1000, " μM")
-confidence_inter = confint(fit; level=0.95)
-
-#---
-plot(xdata, ydata, lab="Data", line=:dash, title="PP1")
-plot!(x-> model_pp1(x, pestim), xdata, lab="Fitted", line=:dot)
-
-# ## Least-square fitting of PLBp
-# Signed Hill function
-hil_s(x, k, n) = sign(x*k) * (abs(x)^n) / (abs(x)^n + abs(k)^n)
-xdata = iso
-ydata = extract(sim, sys.PLBp/sys.PLBtot)
-@. model_plb(x, p) = p[1] * hil_s(x, p[2], p[3]) + 0.1
-p0 = [0.8, 2e-6, 1.0]
-lb = [0.0, 0.0, 0.1]
-fit = curve_fit(model_plb, xdata, ydata, p0; lower=lb)
-pestim = coef(fit)
-
-#---
-confidence_inter = confint(fit; level=0.95)
-
-# It does not fit as good as the previous ones because it's in fact a cubic equation.
-plot(xdata, ydata, lab="Data", line=:dash, title="PLBp")
-plot!(x-> model_plb(x, pestim), xdata, lab="Fitted", line=:dot)
diff --git a/docs/iso.ipynb b/docs/iso.ipynb
deleted file mode 100644
index ac14a82..0000000
--- a/docs/iso.ipynb
+++ /dev/null
@@ -1,215 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Effects of isoproterenol"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "using ModelingToolkit\n",
-    "using DifferentialEquations\n",
-    "using Plots\n",
-    "using CaMKIIModel\n",
-    "Plots.default(fmt=:png)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "sys = build_neonatal_ecc_sys(simplify=true)\n",
-    "tend = 300.0\n",
-    "prob = ODEProblem(sys, [], tend)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## No isoproterenol"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "@unpack Istim = sys\n",
-    "callback = build_stim_callbacks(Istim, tend; period=1)\n",
-    "alg = FBDF()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "sol = solve(prob, alg; callback, abstol=1e-6, reltol=1e-6, maxiters=Int(1e8))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol, idxs=sys.PKACII/sys.RIItot, title=\"PKA activation\", lab=\"PKACII\")\n",
-    "plot!(sol, idxs=sys.PKACI/sys.RItot, title=\"PKA activation\", lab=\"PKACI\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol, idxs=sys.vm, tspan=(295, 300), title=\"Action potential\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol, idxs=[sys.Cai_sub_SR, sys.Cai_sub_SL, sys.Cai_mean], tspan=(299, 300), title=\"Calcium transcient\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol, idxs=sys.CaMKAct, title=\"Active CaMKII\")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## 0.1uM isoproterenol"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "prob2 = remake(prob, p=[sys.ISO => 1E-4])"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "sol2 = solve(prob2, alg; callback, abstol=1e-6, reltol=1e-6, maxiters=Int(1e8))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol2, idxs=sys.PKACII/sys.RIItot, title=\"PKA activation\", lab=\"PKACII\")\n",
-    "plot!(sol2, idxs=sys.PKACI/sys.RItot, title=\"PKA activation\", lab=\"PKACI\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol2, idxs=sys.vm, tspan=(295, 300), title=\"Action potential\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol2, idxs=[sys.Cai_sub_SR, sys.Cai_sub_SL, sys.Cai_mean], tspan=(299, 300), title=\"Calcium transcient\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol2, idxs=sys.CaMKAct, title=\"Active CaMKII\")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Comparison\n",
-    "\n",
-    "Why the opposite effect?"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol, idxs=sys.Cai_mean, tspan=(299, 300), title=\"Calcium transcient\", lab=\"ISO (0uM)\")\n",
-    "plot!(sol2, idxs=sys.Cai_mean, tspan=(299, 300), lab=\"ISO (0.1uM)\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol, idxs=sys.CaMKAct, title=\"Active CaMKII\", lab=\"ISO (0uM)\")\n",
-    "plot!(sol2, idxs=sys.CaMKAct, lab=\"ISO (0.1uM)\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol, idxs=sys.vm, tspan=(297, 300), title=\"Action potential\", lab=\"ISO (0uM)\")\n",
-    "plot!(sol2, idxs=sys.vm, tspan=(297, 300), lab=\"ISO (0.1uM)\")"
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Julia 1.10.5",
-   "language": "julia",
-   "name": "julia-1.10"
-  },
-  "language_info": {
-   "file_extension": ".jl",
-   "mimetype": "application/julia",
-   "name": "julia",
-   "version": "1.10.5"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/docs/iso.jl b/docs/iso.jl
new file mode 100644
index 0000000..e6decf9
--- /dev/null
+++ b/docs/iso.jl
@@ -0,0 +1,59 @@
+# # Effects of isoproterenol
+using ModelingToolkit
+using DifferentialEquations
+using Plots
+using CaMKIIModel
+
+#---
+sys = build_neonatal_ecc_sys(simplify=true)
+tend = 300.0
+prob = ODEProblem(sys, [], tend)
+
+# ## Without isoproterenol
+@unpack Istim = sys
+callback = build_stim_callbacks(Istim, tend; period=1)
+alg = FBDF()
+sol = solve(prob, alg; callback, abstol=1e-6, reltol=1e-6, maxiters=Int(1e8))
+
+#---
+plot(sol, idxs=sys.PKACII/sys.RIItot, title="PKA activation", lab="PKACII")
+plot!(sol, idxs=sys.PKACI/sys.RItot, lab="PKACI")
+
+#---
+plot(sol, idxs=sys.vm, tspan=(295, 300), title="Action potential")
+
+#---
+plot(sol, idxs=[sys.Cai_sub_SR, sys.Cai_sub_SL, sys.Cai_mean], tspan=(299, 300), title="Calcium transcient")
+
+#---
+plot(sol, idxs=sys.CaMKAct, title="Active CaMKII")
+
+# ## 0.1uM isoproterenol
+prob2 = remake(prob, p=[sys.ISO => 1E-4])
+sol2 = solve(prob2, alg; callback, abstol=1e-6, reltol=1e-6, maxiters=Int(1e8))
+
+#---
+plot(sol2, idxs=sys.PKACII/sys.RIItot, title="PKA activation", lab="PKACII")
+plot!(sol2, idxs=sys.PKACI/sys.RItot, lab="PKACI")
+
+#---
+plot(sol2, idxs=sys.vm, tspan=(295, 300), title="Action potential")
+
+#---
+plot(sol2, idxs=[sys.Cai_sub_SR, sys.Cai_sub_SL, sys.Cai_mean], tspan=(299, 300), title="Calcium transcient")
+
+#---
+plot(sol2, idxs=sys.CaMKAct, title="Active CaMKII")
+
+
+# ## Comparison
+plot(sol, idxs=sys.Cai_mean, tspan=(299, 300), title="Calcium transcient", lab="ISO (0uM)")
+plot!(sol2, idxs=sys.Cai_mean, tspan=(299, 300), lab="ISO (0.1uM)")
+
+#---
+plot(sol, idxs=sys.CaMKAct, title="Active CaMKII", lab="ISO (0uM)")
+plot!(sol2, idxs=sys.CaMKAct, lab="ISO (0.1uM)")
+
+#---
+plot(sol, idxs=sys.vm, tspan=(297, 300), title="Action potential", lab="ISO (0uM)")
+plot!(sol2, idxs=sys.vm, tspan=(297, 300), lab="ISO (0.1uM)")
diff --git a/docs/ros.ipynb b/docs/ros.ipynb
deleted file mode 100644
index 607e1d9..0000000
--- a/docs/ros.ipynb
+++ /dev/null
@@ -1,187 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## ROS effects"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "using ModelingToolkit\n",
-    "using DifferentialEquations\n",
-    "using Plots\n",
-    "using CaMKIIModel\n",
-    "Plots.default(fmt=:png)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "sys = build_neonatal_ecc_sys(simplify=true)\n",
-    "tend = 300.0\n",
-    "prob = ODEProblem(sys, [], tend)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## No ROS"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "@unpack Istim = sys\n",
-    "callback = build_stim_callbacks(Istim, tend; period=1)\n",
-    "alg = FBDF()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "sol = solve(prob, alg; callback, abstol=1e-6, reltol=1e-6, maxiters=Int(1e8))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol, idxs=sys.vm, tspan=(298, 300), title=\"Action potential\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol, idxs=[sys.Cai_sub_SR, sys.Cai_sub_SL, sys.Cai_mean], tspan=(299, 300), title=\"Calcium transcient\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol, idxs=sys.CaMKAct, title=\"Active CaMKII\")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## ROS 1uM"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "prob2 = remake(prob, p=[sys.ROS => 1e-4])\n",
-    "sol2 = solve(prob2, alg; callback, abstol=1e-6, reltol=1e-6, maxiters=Int(1e8))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol2, idxs=sys.vm, tspan=(298, 300), title=\"Action potential\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol2, idxs=[sys.Cai_sub_SR, sys.Cai_sub_SL, sys.Cai_mean], tspan=(298, 300), title=\"Calcium transcient\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol2, idxs=sys.CaMKAct, title=\"Active CaMKII\")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Comparisons\n",
-    "\n",
-    "Minimal difference."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol, idxs=sys.CaMKAct, title=\"Active CaMKII\", lab=\"ROS=0\")\n",
-    "plot!(sol2, idxs=sys.CaMKAct, lab=\"ROS=1uM\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol, idxs=sys.vm, title=\"Action potential\", lab=\"ROS=0\", tspan=(298, 300))\n",
-    "plot!(sol2, idxs=sys.vm, lab=\"ROS=1uM\", tspan=(298, 300))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot(sol, idxs=sys.Cai_mean, title=\"Calcium\", lab=\"ROS=0\", tspan=(298, 300))\n",
-    "plot!(sol2, idxs=sys.Cai_mean, lab=\"ROS=1uM\", tspan=(298, 300))"
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Julia 1.10.5",
-   "language": "julia",
-   "name": "julia-1.10"
-  },
-  "language_info": {
-   "file_extension": ".jl",
-   "mimetype": "application/julia",
-   "name": "julia",
-   "version": "1.10.5"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/docs/ros.jl b/docs/ros.jl
new file mode 100644
index 0000000..21d0653
--- /dev/null
+++ b/docs/ros.jl
@@ -0,0 +1,50 @@
+# # ROS effects
+using ModelingToolkit
+using DifferentialEquations
+using Plots
+using CaMKIIModel
+
+#---
+sys = build_neonatal_ecc_sys(simplify=true)
+tend = 300.0
+prob = ODEProblem(sys, [], tend)
+
+# ## No ROS
+@unpack Istim = sys
+callback = build_stim_callbacks(Istim, tend; period=1)
+alg = FBDF()
+sol = solve(prob, alg; callback, abstol=1e-6, reltol=1e-6, maxiters=Int(1e8))
+
+#---
+plot(sol, idxs=sys.vm, tspan=(298, 300), title="Action potential")
+
+#---
+plot(sol, idxs=[sys.Cai_sub_SR, sys.Cai_sub_SL, sys.Cai_mean], tspan=(299, 300), title="Calcium transcient")
+
+#---
+plot(sol, idxs=sys.CaMKAct, title="Active CaMKII")
+
+# ## ROS 1uM
+prob2 = remake(prob, p=[sys.ROS => 1e-4])
+sol2 = solve(prob2, alg; callback, abstol=1e-6, reltol=1e-6, maxiters=Int(1e8))
+
+#---
+plot(sol2, idxs=sys.vm, tspan=(298, 300), title="Action potential")
+
+#---
+plot(sol2, idxs=[sys.Cai_sub_SR, sys.Cai_sub_SL, sys.Cai_mean], tspan=(298, 300), title="Calcium transcient")
+
+#---
+plot(sol2, idxs=sys.CaMKAct, title="Active CaMKII")
+
+# ## Comparisons
+plot(sol, idxs=sys.CaMKAct, title="Active CaMKII", lab="ROS=0")
+plot!(sol2, idxs=sys.CaMKAct, lab="ROS=1uM")
+
+#---
+plot(sol, idxs=sys.vm, title="Action potential", lab="ROS=0", tspan=(298, 300))
+plot!(sol2, idxs=sys.vm, lab="ROS=1uM", tspan=(298, 300))
+
+#---
+plot(sol, idxs=sys.Cai_mean, title="Calcium", lab="ROS=0", tspan=(298, 300))
+plot!(sol2, idxs=sys.Cai_mean, lab="ROS=1uM", tspan=(298, 300))