diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 066c5bf..1df8087 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -13,9 +13,9 @@ concurrency:
 
 env:
   NBCONVERT_JOBS: '2'
-  LITERATE_PROC: '4'
+  LITERATE_PROC: '2'
+  JULIA_NUM_THREADS: '2'
   ALLOWERRORS: 'false'
-  CACHE_NUM: '1'
   JULIA_CONDAPKG_BACKEND: 'Null'
   JULIA_CI: 'true'
   NBCACHE: '.cache'
@@ -23,7 +23,7 @@ env:
 
 jobs:
   CI:
-    runs-on: self-hosted
+    runs-on: ubuntu-latest
     steps:
     - name: Checkout repository
       uses: actions/checkout@v4
@@ -32,8 +32,10 @@ jobs:
       id: setup-python
       with:
         python-version: '3.x'
+    - name: Install uv
+      run: curl -LsSf https://astral.sh/uv/install.sh | sh
     - name: Install Python dependencies
-      run: pip install -r requirements.txt
+      run: uv pip install -r requirements.txt
     - name: Read Julia version
       uses: SebRollen/toml-action@v1.2.0
       id: read_toml
@@ -44,7 +46,7 @@ jobs:
       id: hash
       run: |
         echo "value=${{ hashFiles('Project.toml', 'Manifest.toml', 'src/**') }}" >> "$GITHUB_OUTPUT"
-        echo "ver=${{ runner.os }}-julia-${{ env.CACHE_NUM }}-${{ steps.read_toml.outputs.value }}" >> "$GITHUB_OUTPUT"
+        echo "ver=${{ runner.os }}-julia-${{ steps.read_toml.outputs.value }}" >> "$GITHUB_OUTPUT"
     - name: Cache executed notebooks
       uses: actions/cache@v4
       id: cache-nb
@@ -59,7 +61,7 @@ jobs:
         version: ${{ steps.read_toml.outputs.value }}
     - name: Restore Julia packages
       uses: actions/cache/restore@v4
-      if: ${{ runner.environment == 'github-hosted' }}
+      if: ${{ runner.environment == 'github-hosted'}}
       id: cache-julia
       with:
         path: ~/.julia
@@ -70,7 +72,7 @@ jobs:
       if: ${{ runner.environment == 'self-hosted' || steps.cache-julia.outputs.cache-hit != 'true' }}
       shell: julia --color=yes {0}
       run: |
-        using Pkg, Dates
+        using Pkg
         Pkg.add(["IJulia", "Literate", "PrettyTables", "JSON"])
         Pkg.activate(".")
         Pkg.instantiate()
@@ -83,7 +85,7 @@ jobs:
         key: ${{ steps.cache-julia.outputs.cache-primary-key }}
     - name: Run notebooks
       if: ${{ steps.cache-nb.outputs.cache-hit != 'true' }}
-      run: julia --color=yes -p ${{ env.LITERATE_PROC }} --heap-size-hint=3G run.jl
+      run: julia --project=@. --color=yes -p ${{ env.LITERATE_PROC }} run.jl
     - name: Copy back built notebooks
       run: cp --verbose -rf ${{ env.NBCACHE }}/docs/* docs/
     - name: Build website
diff --git a/.github/workflows/clear-main-cache.yml b/.github/workflows/clear-main-cache.yml
new file mode 100644
index 0000000..05097af
--- /dev/null
+++ b/.github/workflows/clear-main-cache.yml
@@ -0,0 +1,29 @@
+name: Clear main caches
+on:
+  workflow_dispatch:
+
+jobs:
+  cleanup:
+    permissions:
+      actions: write
+    runs-on: ubuntu-latest
+    steps:
+      - name: Cleanup
+        run: |
+          gh extension install actions/gh-actions-cache
+
+          echo "Fetching list of cache key"
+          cacheKeysForPR=$(gh actions-cache list -R $REPO -B $BRANCH -L 100 | cut -f 1 )
+
+          ## Setting this to not fail the workflow while deleting cache keys.
+          set +e
+          echo "Deleting caches..."
+          for cacheKey in $cacheKeysForPR
+          do
+              gh actions-cache delete $cacheKey -R $REPO -B $BRANCH --confirm
+          done
+          echo "Done"
+        env:
+          GH_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          REPO: ${{ github.repository }}
+          BRANCH: ${{ github.event.repository.default_branch }}
diff --git a/.github/workflows/clear-cache.yml b/.github/workflows/clear-pr-cache.yml
similarity index 100%
rename from .github/workflows/clear-cache.yml
rename to .github/workflows/clear-pr-cache.yml
diff --git a/Manifest.toml b/Manifest.toml
index 28e3ac7..67139df 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -2,12 +2,12 @@
 
 julia_version = "1.10.5"
 manifest_format = "2.0"
-project_hash = "c31fe08d86ec107a97ce2f38a0768055cf0cb993"
+project_hash = "dc71554e00e1c1774b67892932418d9d6263c240"
 
 [[deps.ADTypes]]
-git-tree-sha1 = "99a6f5d0ce1c7c6afdb759daa30226f71c54f6b0"
+git-tree-sha1 = "5a5eafb8344b81b8c2237f8a6f6b3602b3f6180e"
 uuid = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
-version = "1.7.1"
+version = "1.8.1"
 weakdeps = ["ChainRulesCore", "EnzymeCore"]
 
     [deps.ADTypes.extensions]
@@ -20,24 +20,28 @@ uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 version = "0.4.5"
 
 [[deps.Accessors]]
-deps = ["CompositionsBase", "ConstructionBase", "Dates", "InverseFunctions", "LinearAlgebra", "MacroTools", "Markdown", "Test"]
-git-tree-sha1 = "f61b15be1d76846c0ce31d3fcfac5380ae53db6a"
+deps = ["CompositionsBase", "ConstructionBase", "InverseFunctions", "LinearAlgebra", "MacroTools", "Markdown"]
+git-tree-sha1 = "b392ede862e506d451fc1616e79aa6f4c673dab8"
 uuid = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
-version = "0.1.37"
+version = "0.1.38"
 
     [deps.Accessors.extensions]
     AccessorsAxisKeysExt = "AxisKeys"
+    AccessorsDatesExt = "Dates"
     AccessorsIntervalSetsExt = "IntervalSets"
     AccessorsStaticArraysExt = "StaticArrays"
     AccessorsStructArraysExt = "StructArrays"
+    AccessorsTestExt = "Test"
     AccessorsUnitfulExt = "Unitful"
 
     [deps.Accessors.weakdeps]
     AxisKeys = "94b1ba4f-4ee9-5380-92f1-94cde586c3c5"
+    Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
     IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
     Requires = "ae029012-a4dd-5104-9daa-d747884805df"
     StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
     StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
+    Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
     Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [[deps.Adapt]]
@@ -130,9 +134,9 @@ version = "0.1.6"
 
 [[deps.BlockArrays]]
 deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra"]
-git-tree-sha1 = "5c0ffe1dff8cb7112de075f1b1cb32191675fcba"
+git-tree-sha1 = "d434647f798823bcae510aee0bc0401927f64391"
 uuid = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
-version = "1.1.0"
+version = "1.1.1"
 
     [deps.BlockArrays.extensions]
     BlockArraysBandedMatricesExt = "BandedMatrices"
@@ -166,9 +170,9 @@ version = "1.18.0+2"
 
 [[deps.ChainRulesCore]]
 deps = ["Compat", "LinearAlgebra"]
-git-tree-sha1 = "71acdbf594aab5bbb2cec89b208c41b4c411e49f"
+git-tree-sha1 = "3e4b134270b372f2ed4d4d0e936aabaefc1802bc"
 uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-version = "1.24.0"
+version = "1.25.0"
 weakdeps = ["SparseArrays"]
 
     [deps.ChainRulesCore.extensions]
@@ -342,9 +346,9 @@ version = "1.9.1"
 
 [[deps.DiffEqBase]]
 deps = ["ArrayInterface", "ConcreteStructs", "DataStructures", "DocStringExtensions", "EnumX", "EnzymeCore", "FastBroadcast", "FastClosures", "ForwardDiff", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "Logging", "Markdown", "MuladdMacro", "Parameters", "PreallocationTools", "PrecompileTools", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SciMLStructures", "Setfield", "Static", "StaticArraysCore", "Statistics", "Tricks", "TruncatedStacktraces"]
-git-tree-sha1 = "d7d43a1cc11dc4e4e5378816ae720fccd75a77c8"
+git-tree-sha1 = "fa7d580038451a8df4434ef5b079ac9b2d486194"
 uuid = "2b5f629d-d688-5b77-993f-72d75c75574e"
-version = "6.155.0"
+version = "6.155.1"
 
     [deps.DiffEqBase.extensions]
     DiffEqBaseCUDAExt = "CUDA"
@@ -904,9 +908,9 @@ version = "3.0.0+1"
 
 [[deps.LLVMOpenMP_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "e16271d212accd09d52ee0ae98956b8a05c4b626"
+git-tree-sha1 = "78211fb6cbc872f77cad3fc0b6cf647d923f4929"
 uuid = "1d63c593-3942-5779-bab2-d838dc0a180e"
-version = "17.0.6+0"
+version = "18.1.7+0"
 
 [[deps.LZO_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -1189,9 +1193,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", "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 = "bc5164aa274590945239c84b198fd5ce475a4772"
+git-tree-sha1 = "1f63c7127bbabf22c58edc85dac72afa28fce375"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
-version = "9.39.1"
+version = "9.40.0"
 
     [deps.ModelingToolkit.extensions]
     MTKBifurcationKitExt = "BifurcationKit"
@@ -1671,6 +1675,12 @@ version = "0.5.6"
 deps = ["Unicode"]
 uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
+[[deps.ProgressLogging]]
+deps = ["Logging", "SHA", "UUIDs"]
+git-tree-sha1 = "80d919dee55b9c50e8d9e2da5eeafff3fe58b539"
+uuid = "33c8b6b6-d38a-422a-b730-caa89a2f386c"
+version = "0.1.4"
+
 [[deps.PtrArrays]]
 git-tree-sha1 = "77a42d78b6a92df47ab37e177b2deac405e1c88f"
 uuid = "43287f4e-b6f4-7ad1-bb20-aadabca52c3d"
@@ -1702,9 +1712,9 @@ version = "6.7.1+1"
 
 [[deps.QuadGK]]
 deps = ["DataStructures", "LinearAlgebra"]
-git-tree-sha1 = "1d587203cf851a51bf1ea31ad7ff89eff8d625ea"
+git-tree-sha1 = "cda3b045cf9ef07a08ad46731f5a3165e56cf3da"
 uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
-version = "2.11.0"
+version = "2.11.1"
 
     [deps.QuadGK.extensions]
     QuadGKEnzymeExt = "Enzyme"
@@ -1801,15 +1811,15 @@ version = "1.1.1"
 
 [[deps.Rmath]]
 deps = ["Random", "Rmath_jll"]
-git-tree-sha1 = "f65dcb5fa46aee0cf9ed6274ccbd597adc49aa7b"
+git-tree-sha1 = "852bd0f55565a9e973fcfee83a84413270224dc4"
 uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
-version = "0.7.1"
+version = "0.8.0"
 
 [[deps.Rmath_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "e60724fd3beea548353984dc61c943ecddb0e29a"
+git-tree-sha1 = "58cdd8fb2201a6267e1db87ff148dd6c1dbd8ad8"
 uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
-version = "0.4.3+0"
+version = "0.5.1+0"
 
 [[deps.RuntimeGeneratedFunctions]]
 deps = ["ExprTools", "SHA", "Serialization"]
@@ -1834,9 +1844,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 = "952e09c65264fe2a8155dd6aad60c188a0b83f69"
+git-tree-sha1 = "c96f81c3e98d5e2f0848fb42aba4383d772c3bb7"
 uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
-version = "2.53.0"
+version = "2.53.1"
 
     [deps.SciMLBase.extensions]
     SciMLBaseChainRulesCoreExt = "ChainRulesCore"
@@ -1906,9 +1916,9 @@ version = "1.1.0"
 
 [[deps.SimpleNonlinearSolve]]
 deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "DiffEqBase", "DiffResults", "DifferentiationInterface", "FastClosures", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "MaybeInplace", "PrecompileTools", "Reexport", "SciMLBase", "Setfield", "StaticArraysCore"]
-git-tree-sha1 = "4d7a7c177bcb4c6dc465f09db91bfdb28c578919"
+git-tree-sha1 = "536c0ee0b4b766ddee24220c6bb60932df4e2c39"
 uuid = "727e6d20-b764-4bd8-a329-72de5adea6c7"
-version = "1.12.0"
+version = "1.12.1"
 
     [deps.SimpleNonlinearSolve.extensions]
     SimpleNonlinearSolveChainRulesCoreExt = "ChainRulesCore"
@@ -2041,9 +2051,9 @@ version = "0.34.3"
 
 [[deps.StatsFuns]]
 deps = ["HypergeometricFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
-git-tree-sha1 = "cef0472124fab0695b58ca35a77c6fb942fdab8a"
+git-tree-sha1 = "b423576adc27097764a90e163157bcfc9acf0f46"
 uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
-version = "1.3.1"
+version = "1.3.2"
 weakdeps = ["ChainRulesCore", "InverseFunctions"]
 
     [deps.StatsFuns.extensions]
@@ -2079,9 +2089,9 @@ version = "0.2.2"
 
 [[deps.SymbolicUtils]]
 deps = ["AbstractTrees", "ArrayInterface", "Bijections", "ChainRulesCore", "Combinatorics", "ConstructionBase", "DataStructures", "DocStringExtensions", "DynamicPolynomials", "IfElse", "LinearAlgebra", "MultivariatePolynomials", "NaNMath", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArrays", "SymbolicIndexingInterface", "TermInterface", "TimerOutputs", "Unityper"]
-git-tree-sha1 = "9d983078d9e99421fcca44c373e4304b8421fdde"
+git-tree-sha1 = "3927e02dc7648a45ec6aa592bcd8374094a44740"
 uuid = "d1185830-fcd6-423d-90d6-eec64667417b"
-version = "3.6.0"
+version = "3.7.1"
 
     [deps.SymbolicUtils.extensions]
     SymbolicUtilsLabelledArraysExt = "LabelledArrays"
@@ -2093,9 +2103,9 @@ version = "3.6.0"
 
 [[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 = "2226d810512c678d2ec9c2a9b2e227c2ebc43573"
+git-tree-sha1 = "8b48697e7fec6d4b7c4a9fe892857a5ed2bae7e8"
 uuid = "0c5d862f-8b57-4792-8d23-62f2024744c7"
-version = "6.11.0"
+version = "6.12.0"
 
     [deps.Symbolics.extensions]
     SymbolicsForwardDiffExt = "ForwardDiff"
diff --git a/Project.toml b/Project.toml
index ef40b37..587f921 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "CaMKIIModel"
 uuid = "268d5184-1524-48b5-bbe3-a3af95dac460"
 authors = ["Wen-Wei Tseng <sosiristseng@gmail.com>"]
-version = "0.2.0"
+version = "0.3.0"
 
 [deps]
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
@@ -9,3 +9,4 @@ ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+ProgressLogging = "33c8b6b6-d38a-422a-b730-caa89a2f386c"
diff --git a/docs/camkii_calwave.jl b/docs/camkii_calwave.jl
index 717c938..b794992 100644
--- a/docs/camkii_calwave.jl
+++ b/docs/camkii_calwave.jl
@@ -10,7 +10,7 @@ using CaMKIIModel: get_camkii_sys, μM, nM, second
 @parameters ROS=0μM period=1/3 ca_r=100nM ca_rise=550nM tstart=200.0second tend=300.0second
 @independent_variables t
 @variables Ca(t)
-sys = get_camkii_sys(;Ca, ROS)
+sys = get_camkii_sys(Ca;ROS)
 
 # Periodic assymetric calcium pulses
 function ca_wave(t;
diff --git a/docs/camkii_ros.jl b/docs/camkii_ros.jl
index 6ffc1ca..529bda7 100644
--- a/docs/camkii_ros.jl
+++ b/docs/camkii_ros.jl
@@ -29,7 +29,7 @@ end
 @independent_variables t
 @variables Ca(t)
 casys = ca_decay_sys()
-sys = get_camkii_sys(;Ca, ROS)
+sys = get_camkii_sys(Ca; ROS)
 sys = extend(casys, sys)
 sys = structural_simplify(sys)
 unknowns(sys)
diff --git a/run.jl b/run.jl
index 31f94c5..16fb7d6 100644
--- a/run.jl
+++ b/run.jl
@@ -6,7 +6,6 @@ using IJulia
 @everywhere begin
     ENV["GKSwstype"] = "100"
     using Literate, Pkg, JSON
-    Pkg.activate(Base.current_project())
 end
 
 # Strip SVG output from a Jupyter notebook
@@ -100,9 +99,7 @@ function main(;
     ts_lit = pmap(litnbs; on_error=ex -> NaN) do nb
         @elapsed run_literate(nb, cachedir; rmsvg)
     end
-
-    # Remove worker processes to release some memory
-    rmprocs(workers())
+    rmprocs(workers()) # Remove worker processes to release some memory
 
     # Debug notebooks one by one if there are errors
     for (nb, t) in zip(litnbs, ts_lit)
diff --git a/src/CaMKIIModel.jl b/src/CaMKIIModel.jl
index 04a09ec..ca25b81 100644
--- a/src/CaMKIIModel.jl
+++ b/src/CaMKIIModel.jl
@@ -3,8 +3,10 @@ module CaMKIIModel
 using ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
 using NaNMath
+using OrdinaryDiffEq
+using DiffEqCallbacks
 
-export get_camkii_sys, get_bar_sys, build_neonatal_ecc_sys
+export get_camkii_sys, get_bar_sys, build_neonatal_ecc_sys, build_stim_callbacks
 
 include("utils.jl")
 include("camkii_ros.jl")
diff --git a/src/camkii_ros.jl b/src/camkii_ros.jl
index 5481b3e..e50d396 100644
--- a/src/camkii_ros.jl
+++ b/src/camkii_ros.jl
@@ -1,6 +1,6 @@
 # CaMKII system with ROS activation
-function get_camkii_sys(;
-    Ca=0μM, ROS=0μM,
+function get_camkii_sys(Ca=0μM;
+    ROS=0μM,
     binding_To_PCaMK=0.1,
     decay_CaM=3,
     phospho_rate=1Hz,
@@ -92,7 +92,6 @@ function get_camkii_sys(;
         CAM_T ~ CaM0 + Ca2CaM_C + Ca2CaM_N + Ca4CaM + CaM0_CaMK + Ca2CaM_C_CaMK + Ca2CaM_N_CaMK + Ca4CaM_CaMK + CaM0_CaMKP + Ca2CaM_C_CaMKP + Ca2CaM_N_CaMKP + Ca4CaM_CaMKP + Ca4CaM_CaMKOX + Ca4CaM_CaMKPOX,
     ]
 
-    Num(0)
     rates = merge(Dict(sts .=> Num(0)), Dict(conservedvars .=> Num(0)))
 
     # Observables
diff --git a/src/capde.jl b/src/capde.jl
index ebbe6ce..d4b45bf 100644
--- a/src/capde.jl
+++ b/src/capde.jl
@@ -1,28 +1,13 @@
-# Calcium diffusion between sarcolemma (SL) and sarcoplasmic reticulum (SR)
-"Calcium buffering factor by troponin and calmodulin"
-function beta_cai(Ca;
-    TnI_PKAp=0,
-    TrpnTotal=35μM,
-    CmdnTotal=50μM,
-    KmCmdn=2.38μM,
-    KmTrpn=0.5μM,
-    fracTnIp0=0.062698 # Baseline effect
-)
-    fPKA_TnI = 1.61 - 0.61 * (1 - TnI_PKAp) / (1 - fracTnIp0) # Max effect +61%
-    KmTrpnNew = KmTrpn / fPKA_TnI
-    return inv(1 + TrpnTotal * KmTrpnNew / (Ca + KmTrpnNew)^2 + CmdnTotal * KmCmdn / (Ca + KmCmdn)^2)
-end
-
 "Calcium diffusion between sarcolemma (SL) and sarcoplasmic reticulum (SR)"
 function get_ca_pde_sys(;
-    Cai_sub_SR_default=0.2556μM,
-    Cai_sub_SL_default=0.25151μM,
+    Cai_default=0.2556μM,
     dx=0.1μm,
     rSR_true=6μm,
     rSL_true=10.5μm,
     TnI_PKAp=0,
     JCa_SR=0,
     JCa_SL=0,
+    V_sub_SL=0.137pL,
     name=:capdesys
 )
     rSR = rSR_true + 0.5 * dx
@@ -31,29 +16,33 @@ function get_ca_pde_sys(;
     m = length(j)
     @variables Cai(t)[1:m] Cai_mean(t) Cai_sub_SR(t) Cai_sub_SL(t)
     @parameters begin
-        Dca = 7 * (μm^2 // ms)
-        V_sub_SR = 4 // 3 * pi * ((rSR_true + dx)^3 - (rSR_true)^3)
-        V_sub_SL = 4 // 3 * pi * (rSL_true^3 - (rSL_true - dx)^3)
+        Dca = 7μm^2 / ms
         TrpnTotal = 35μM
         CmdnTotal = 50μM
         KmCmdn = 2.38μM
         KmTrpn = 0.5μM
         fracTnIp0 = 0.062698 # Baseline effect
     end
+
+    fPKA_TnI = 1.61 - 0.61 * (1 - TnI_PKAp) / (1 - fracTnIp0) # Max effect +61%
+    KmTrpnNew = KmTrpn / fPKA_TnI
+
+    "Calcium buffering factor by troponin and calmodulin"
+    beta_cai(Ca) = inv(1 + TrpnTotal * KmTrpnNew / (Ca + KmTrpnNew)^2 + CmdnTotal * KmCmdn / (Ca + KmCmdn)^2)
+
     eqs = [
         Cai_mean ~ sum(collect(Cai)) / m,
         Cai_sub_SR ~ Cai[1],
         Cai_sub_SL ~ Cai[m],
-        D(Cai[1]) ~ (Dca / (j[1] * dx^2) * ((1 + j[1]) * Cai[2] - 2 * j[1] * Cai[1] + (j[1] - 1) * Cai[1]) + JCa_SR / V_sub_SR) * beta_cai(Cai[1]; TnI_PKAp, TrpnTotal, CmdnTotal, KmCmdn, KmTrpn, fracTnIp0),
-        D(Cai[m]) ~ (Dca / (j[m] * dx^2) * ((1 + j[m]) * Cai[m] - 2 * j[m] * Cai[m] + (j[m] - 1) * Cai[m-1]) + JCa_SL / V_sub_SL) * beta_cai(Cai[m]; TnI_PKAp, TrpnTotal, CmdnTotal, KmCmdn, KmTrpn, fracTnIp0),
+        D(Cai[1]) ~ (Dca / (j[1] * dx^2) * ((1 + j[1]) * Cai[2] - 2 * j[1] * Cai[1] + (j[1] - 1) * Cai[1]) + JCa_SR) * beta_cai(Cai[1]),
+        D(Cai[m]) ~ (Dca / (j[m] * dx^2) * ((1 + j[m]) * Cai[m] - 2 * j[m] * Cai[m] + (j[m] - 1) * Cai[m-1]) + JCa_SL / V_sub_SL) * beta_cai(Cai[m]),
     ]
 
-    defaults = [Cai[1] => Cai_sub_SR_default, Cai[m] => Cai_sub_SL_default]
+    defaults = [Cai[1] => Cai_default, Cai[m] => Cai_default]
 
     for i in 2:m-1
-        eq = D(Cai[i]) ~ (Dca / (j[i] * dx^2) * ((1 + j[i]) * Cai[i+1] - 2 * j[i] * Cai[i] + (j[i] - 1) * Cai[i-1])) * beta_cai(Cai[i]; TnI_PKAp, TrpnTotal, CmdnTotal, KmCmdn, KmTrpn, fracTnIp0)
-        push!(eqs, eq)
-        push!(defaults, Cai[i] => (Cai_sub_SR_default + Cai_sub_SL_default) / 2)
+        push!(eqs, D(Cai[i]) ~ (Dca / (j[i] * dx^2) * ((1 + j[i]) * Cai[i+1] - 2 * j[i] * Cai[i] + (j[i] - 1) * Cai[i-1])) * beta_cai(Cai[i]))
+        push!(defaults, Cai[i] => Cai_default)
     end
 
     return ODESystem(eqs, t; name, defaults)
diff --git a/src/ecc_neonatal.jl b/src/ecc_neonatal.jl
index 37c3d32..e3cf521 100644
--- a/src/ecc_neonatal.jl
+++ b/src/ecc_neonatal.jl
@@ -5,7 +5,7 @@
 # PLEASE MENTION THE FOLLOWING REFERENCE WHEN USING THIS CODE OR PART OF IT:
 # Korhonen et al. "Model of excitation-contraction coupling of rat neonatal
 # ventricular myocytes" Biophys J. 2009, Feb; 96(3):1189-1209
-function get_nak_sys(vm, Nai, Nao, Ko; name=:naksys)
+function get_nak_sys(na_i, na_o, k_o, vm; name=:naksys)
     @parameters begin
         INaKmax = 2.7μAμF
         KmNaiNaK = 18.6mM
@@ -14,67 +14,71 @@ function get_nak_sys(vm, Nai, Nao, Ko; name=:naksys)
     end
 
     @variables INaK(t)
-    sigma = 1 / 7 * expm1(Nao / 67.300mM)
+    sigma = 1 / 7 * expm1(na_o / 67.3mM)
     fNaK = inv(1 + 0.1245 * exp(-0.1vm * iVT) + 0.0365 * sigma * exp(-vm * iVT))
-    fKo = hil(Ko, KmKoNaK)
-    fNai = hil(Nai, KmNaiNaK, nNaK)
+    fKo = hil(k_o, KmKoNaK)
+    fNai = hil(na_i, KmNaiNaK, nNaK)
     return ODESystem([INaK ~ INaKmax * fNaK * fKo * fNai], t; name)
 end
 
 function build_neonatal_ecc_sys(;
     rSR_true=6μm,
     rSL_true=10.5μm,
+    dx=0.1μm,
     name=:neonataleccsys,
     simplify=true,
 )
     @parameters begin
-        Ca_o = 1796μM
-        Na_o = 154578μM
-        K_o = 5366μM
-        Mg_i = 1000μM
+        ca_o = 1796μM
+        na_o = 154578μM
+        k_o = 5366μM
+        mg_i = 1000μM
         ROS = 0μM
         ISO = 0μM
         ATP = 5mM
-        Cm = 1μF // cm^2
-        Acap = 4π * rSL_true^2
-        Vmyo = 4 // 3 * π * (rSL_true^3 - rSR_true^3)
-        ACAP_VMYO_F = Acap * Cm / Faraday / Vmyo
         Istim = 0
+        # cell geometry
+        Cm = 1μF / cm^2
+        Acap = 4π * rSL_true^2
+        Vmyo = 4 / 3 * π * (rSL_true^3 - rSR_true^3) # 3.944 pL
+        ACAP_F = Acap * Cm / Faraday
+        V_sub_SR = 4 / 3 * pi * ((rSR_true + dx)^3 - (rSR_true)^3) # 0.046 pL
+        V_sub_SL = 4 / 3 * pi * (rSL_true^3 - (rSL_true - dx)^3)   # 0.137 pL
     end
 
     @variables begin
-        Na_i(t) = 13838.37602μM
-        K_i(t) = 150952.75035μM
+        na_i(t) = 13838.37602μM
+        k_i(t) = 150952.75035μM
         vm(t) = -68.79268mV
         JCa_SL(t)
         JCa_SR(t)
     end
 
-    barsys = get_bar_sys(; ATP, ISO)
+    barsys = get_bar_sys(ATP, ISO)
     @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)
+    capdesys = get_ca_pde_sys(; JCa_SR, JCa_SL, TnI_PKAp, rSR_true, rSL_true, dx, V_sub_SL)
     @unpack Cai_sub_SL, Cai_sub_SR, Cai_mean = capdesys
-    camkiisys = get_camkii_sys(; ROS, Ca=Cai_mean)
-    icasys = get_ica_sys(Na_i, Cai_sub_SL, Na_o, Ca_o, vm; LCCb_PKAp)
+    camkiisys = get_camkii_sys(Cai_mean; ROS)
+    icasys = get_ica_sys(na_i, Cai_sub_SL, na_o, ca_o, vm; Acap, Cm, LCCb_PKAp)
     @unpack INaCa, ICaL, ICaT, ICab = icasys
-    inasys = get_ina_sys(Na_i, Na_o, vm)
+    inasys = get_ina_sys(na_i, na_o, vm)
     @unpack INa, INab = inasys
-    iksys = get_ik_sys(K_i, K_o, Na_i, Na_o, vm; IKUR_PKAp)
+    iksys = get_ik_sys(k_i, k_o, na_i, na_o, vm; IKUR_PKAp)
     @unpack IK1, Ito, IKs, IKr, IfNa, IfK, If = iksys
-    sersys = get_ser_sys(Cai_sub_SR; fracPLBp)
-    naksys = get_nak_sys(vm, Na_i, Na_o, K_o)
+    sersys = get_ser_sys(Cai_sub_SR; fracPLBp, V_sub_SR)
+    naksys = get_nak_sys(na_i, na_o, k_o, vm)
     @unpack INaK = naksys
 
-    sys = ODESystem([
-        D(vm) ~ -(INab + INaCa + ICaL + ICaT + If + Ito + IK1 + IKs + IKr + INa + INaK + ICab + Istim), # Current alreadt normalized by capacitance
-        D(Na_i) ~ -(IfNa + INab + INa + 3 * INaCa + 3 * INaK) * ACAP_VMYO_F,
-        D(K_i) ~ -(IfK + Ito + IK1 + IKs + IKr + Istim - 2 * INaK) * ACAP_VMYO_F,
-    ], t; name)
+    eqs = [
+        D(vm) ~ -(INab + INaCa + ICaL + ICaT + If + Ito + IK1 + IKs + IKr + INa + INaK + ICab + Istim), # Current  normalized by capacitance
+        D(na_i) ~ -(IfNa + INab + INa + 3 * INaCa + 3 * INaK) * ACAP_F / Vmyo,
+        D(k_i) ~ -(IfK + Ito + IK1 + IKs + IKr + Istim - 2 * INaK) * ACAP_F / Vmyo,
+    ]
 
+    sys = ODESystem(eqs, t; name)
     for s2 in (barsys, capdesys, camkiisys, icasys, inasys, iksys, sersys, naksys)
-        sys = extend(s2, sys; name)
+        sys = extend(sys, s2; name)
     end
-
     if simplify
         sys = structural_simplify(sys)
     end
diff --git a/src/er.jl b/src/er.jl
index 40b1b46..b08c60c 100644
--- a/src/er.jl
+++ b/src/er.jl
@@ -1,24 +1,20 @@
 "Sarcoplasmic reticulum system"
-function get_ser_sys(Cai; fracPLB_CKp=0, fracPLBp=0, RyR_CKp=0, name=:sersys)
+function get_ser_sys(Cai_sub_SR; fracPLB_CKp=0, fracPLBp=0, RyR_CKp=0.2, V_sub_SR=0.046pL, name=:sersys)
     @parameters begin
-        VSR = 0.043 * 1.5 * 1.4pL
+        VSR = 0.093pL
         VNSR = 0.9 * VSR
         VJSR = VSR - VNSR
         # RyR
-        nRyR = 4
-        nu1RyR = 0.01 / ms
-        kaposRyR = 1 / ms
-        kanegRyR = 0.16 / ms
+        kRyR = 10Hz
+        kaposRyR = 1000Hz
+        kanegRyR = 160Hz
         # SERCA
-        VmaxfSR = 0.9996 * (μM / ms)
-        VmaxrSR = VmaxfSR
+        VmaxSR = 0.9996 * μM / ms
         KmfSR = 0.5μM
         KmrSR = 7000 * KmfSR
-        HfSR = 2
-        HrSR = 1 * HfSR
         kSRleak = 5e-6 / ms
         fracPKA_PLBo = 1 - 0.079755
-        ktrCaSR = 1 / 200ms
+        ktrCaSR = inv(200ms)
         csqntot = 24750μM
         Kmcsqn = 800μM
     end
@@ -34,27 +30,31 @@ function get_ser_sys(Cai; fracPLB_CKp=0, fracPLBp=0, RyR_CKp=0, name=:sersys)
         Jtr(t)
         betaSR(t)
         JCa_SR(t)
+        KmRyR(t)
+        dPO1RyR(t)
     end
 
-    KmRyR = (1.35 * 2.6 * expit(-(Cai - 530μM) / 200μM) + 1.5 - 0.9 - 0.3 - 0.05) * μM
     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 smaller value (resulting in max reduction of Kmf)
+    # Select the smaller value (resulting in max reduction of Kmf)
     Kmfp = 2 * KmfSR * min(fCKII_PLB, fPKA_PLB)  #fCKII_PLB
-    fSR = NaNMath.pow(Cai / Kmfp, HfSR)
-    rSR = NaNMath.pow(CaNSR / KmrSR, HrSR)
-    kleak = (1 / 2 + 5 * RyR_CKp / 2) * kSRleak
+    fSR = (Cai_sub_SR / Kmfp)^2
+    rSR = (CaNSR / KmrSR)^2
+    kleak = (1 + 5 * RyR_CKp) * kSRleak / 2
 
-    return ODESystem([
+    eqs = [
         1 ~ PO1RyR + PC1RyR,
-        Jrel ~ nu1RyR * PO1RyR * (CaJSR - Cai),
-        D(PO1RyR) ~ kaposRyR * hil(Cai, KmRyR, nRyR) * PC1RyR - kanegRyR * PO1RyR,
-        Jup ~ (VmaxfSR * fSR - VmaxrSR * rSR) / (1 + fSR + rSR),
-        Jleak ~ kleak * (CaNSR - Cai),
+        Jrel ~ kRyR * PO1RyR * (CaJSR - Cai_sub_SR),
+        KmRyR ~ (3.51 / (1 + exp((CaJSR - 530μM) / 200μM)) + 0.25) * μM,
+        dPO1RyR ~ kaposRyR * hil(Cai_sub_SR, KmRyR, 4) * PC1RyR - kanegRyR * PO1RyR,
+        D(PO1RyR) ~ dPO1RyR,
+        Jup ~ VmaxSR * (fSR - rSR) / (1 + fSR + rSR),
+        Jleak ~ kleak * (CaNSR - Cai_sub_SR),
         Jtr ~ ktrCaSR * (CaNSR - CaJSR),
         betaSR ~ inv(1 + csqntot * Kmcsqn / (CaJSR + Kmcsqn)^2),
         JCa_SR ~ Jleak - Jup + Jrel,
-        D(CaJSR) ~ betaSR * (-Jrel + Jtr) / VJSR,
-        D(CaNSR) ~ (Jup - Jleak - Jtr) / VNSR,
-    ], t; name)
+        D(CaJSR) ~ betaSR * (-Jrel * V_sub_SR + Jtr * VNSR)/VJSR,
+        D(CaNSR) ~ (Jup - Jleak) * V_sub_SR / VNSR - Jtr,
+    ]
+    return ODESystem(eqs, t; name)
 end
diff --git a/src/ica.jl b/src/ica.jl
index a07d790..4236b84 100644
--- a/src/ica.jl
+++ b/src/ica.jl
@@ -1,5 +1,5 @@
-# Plama membrane calcium currents
-function get_ica_sys(nai, cai, nao, cao, vm; Acap=4π * (10μm)^2, Cm=1μF // cm^2, LCCb_PKAp=0, name=:icasys)
+"Plama membrane calcium currents"
+function get_ica_sys(nai, cai, nao, cao, vm; Acap=4π * (10μm)^2, Cm=1μF / cm^2, LCCb_PKAp=0, name=:icasys)
     @parameters begin
         ICa_scale0 = 0.95 # or 5.25
         fracLCCbp0 = 0.250657 # Derived quantity - (LCCbp(baseline)/LCCbtot)
@@ -8,7 +8,7 @@ function get_ica_sys(nai, cai, nao, cao, vm; Acap=4π * (10μm)^2, Cm=1μF // cm
         kNaCa = 2.268e-016μAμF / μM^4
         dNaCa = 1e-16 / μM^4
         gNaCa = 0.5
-        GCaL = 6.3e-5 * ((metre//10)^3 // ms // Farad)
+        GCaL = 6.3e-5 * ((0.1metre)^3 / ms / Farad)
         taufca = 10ms
         gCaT = 0.2mSμF
         gCab = 0.0008mSμF
@@ -38,7 +38,7 @@ function get_ica_sys(nai, cai, nao, cao, vm; Acap=4π * (10μm)^2, Cm=1μF // cm
         ICab(t)
     end
 
-    # Calcium flux scaled by LCC (phosphorylated@beta subunit)
+    # Calcium flux scaled by LCC (phosphorylated beta subunit)
     a_favail = (1.56 - 1) / (fracLCCbpISO / fracLCCbp0 - 1)         # fracLCCbp ISO (x1.56 o.1 ISO)
     favail = (1 - a_favail) + a_favail * (LCCb_PKAp / fracLCCbp0)   # Test (max x2.52 100# phosph)
 
@@ -46,7 +46,7 @@ function get_ica_sys(nai, cai, nao, cao, vm; Acap=4π * (10μm)^2, Cm=1μF // cm
     a = nai^3 * cao
     b = nao^3 * cai * fNaCa
     # L-type calcium channel (LCC)
-    V = vm / mV # Convert voltage to mV
+    V = vm * inv(mV) # Convert voltage to mV
     alphad = 1.4 * expit((V + 35) / 13) + 0.25
     betad = 1.4 * expit(-(V + 5) / 5)
     gammad = expit((V - 50) / 20)
@@ -58,7 +58,7 @@ function get_ica_sys(nai, cai, nao, cao, vm; Acap=4π * (10μm)^2, Cm=1μF // cm
     return ODESystem([
             ICa_scale ~ ICa_scale0 * favail,
             E_Ca ~ nernst(cao, cai, 2),
-            JCa_SL ~ (2 * INaCa - ICaL - ICaT - ICab) * (Acap * Cm // Faraday),
+            JCa_SL ~ (2 * INaCa - ICaL - ICaT - ICab) * (Acap * Cm / 2 / Faraday),
             INaCa ~ ICa_scale * kNaCa * (exp(iVT * gNaCa * vm) * a - exp(iVT * (gNaCa - 1) * vm) * b) / (1 + dNaCa * (a + b)),
             ICaL ~ ICa_scale * i_d * i_f * i_fca * ghkVm(GCaL, vm, cai, 0.341 * cao, 2),
             dinf ~ expit((V + 11.1) / 7.2),
diff --git a/src/ik.jl b/src/ik.jl
index f4347b1..62356eb 100644
--- a/src/ik.jl
+++ b/src/ik.jl
@@ -1,8 +1,8 @@
 "Potassium currents"
-function get_ik_sys(K_i, K_o, Na_i, Na_o, vm; IKUR_PKAp=0, name=:iksys)
+function get_ik_sys(k_i, k_o, na_i, na_o, vm; IKUR_PKAp=0, name=:iksys)
     @parameters begin
         # IK1: time-independent
-        gK1 = 0.0515mSμF * hil(K_o, 210μM)
+        gK1 = 0.0515mSμF * hil(k_o, 210μM)
         # Ito: Where does this come from? Perhaps here: https://modeldb.science/262081
         gt = 0.1mSμF
         f_is = 0.706
@@ -49,7 +49,7 @@ function get_ik_sys(K_i, K_o, Na_i, Na_o, vm; IKUR_PKAp=0, name=:iksys)
         i_y(t) = 0.07192
     end
 
-    V = vm * Volt / mV # Convert voltage to mV
+    V = vm * inv(mV)  # Convert voltage to mV
     vk1 = vm - E_K - 6.1373mV
 
     # PKA-dependent phosphoregulation of Ik,slow1 (increases Gkur1)
@@ -61,12 +61,12 @@ function get_ik_sys(K_i, K_o, Na_i, Na_o, vm; IKUR_PKAp=0, name=:iksys)
     betan = 0.0000953333 * exp(-0.038 * (V + 26.5))
 
     # IKr
-    alphaa0 = 0.022348 / ms * exp(0.01176 * V)
-    betaa0 = 0.047002 / ms * exp(-0.0631 * V)
-    alphaa1 = 0.013733 / ms * exp(0.038198 * V)
-    betaa1 = 0.0000689 / ms * exp(-0.04178 * V)
-    alphai_mERG = 0.090821 / ms * exp(0.023391 * V)
-    betai_mERG = 0.006497 / ms * exp(-0.03268 * V)
+    alphaa0 = 0.022348kHz * exp(0.01176 * V)
+    betaa0 = 0.047002kHz * exp(-0.0631 * V)
+    alphaa1 = 0.013733kHz * exp(0.038198 * V)
+    betaa1 = 0.0000689kHz * exp(-0.04178 * V)
+    alphai_mERG = 0.090821kHz * exp(0.023391 * V)
+    betai_mERG = 0.006497kHz * exp(-0.03268 * V)
 
     rc0c1 = alphaa0 * CK0 - betaa0 * i_CK1
     rc1c2 = kf * i_CK1 - kb * i_CK2
@@ -76,34 +76,34 @@ function get_ik_sys(K_i, K_o, Na_i, Na_o, vm; IKUR_PKAp=0, name=:iksys)
     fK = 1 - fNa
 
     return ODESystem([
-        E_Na ~ nernst(Na_o, Na_i),
-        E_K ~ nernst(K_o, K_i),
-        IK1 ~ gK1 * vk1 / (0.1653 + exp(0.0319 / mV * vk1)),
-        sinf ~ expit((vm + 31.97156mV) / -4.64291mV),
-        rinf ~ expit((vm - 3.55716mV) / 14.61299mV),
-        slowinf ~ sinf,
-        taur ~ inv(45.16 * exp(0.03577 / mV * (vm + 50mV)) + 98.9 * exp(-0.1 / mV * (vm + 38mV))),
-        taus ~ (0.35 * exp(-(((vm + 70mV) / 15mV)^2)) + 0.035) - 26.9ms,
-        tausslow ~ (3.7 * exp(-(((vm + 70mV) / 30mV)^2)) + 0.035) + 37.4ms,
-        Ito ~ gt * i_r * (f_is * i_s + (1 - f_is) * i_sslow) * (vm - E_K),
-        D(i_r) ~ (rinf - i_r)/taur,
-        D(i_s) ~ (sinf - i_s)/taus,
-        D(i_sslow) ~ (slowinf - i_sslow)/tausslow,
-        IKs ~ GKs * i_nKs^2 * (vm - E_K) * fracIKuravail * 2,
-        D(i_nKs) ~ (alphan / (alphan + betan) - i_nKs)/nKstau,
-        E_Kr ~ nernst(0.98 * K_o + 0.02 * Na_o, 0.98 * K_i + 0.02 * Na_i),
-        IKr ~ GKr * i_OK * (vm - E_Kr),
-        1 ~ CK0 + i_CK1 + i_CK2 + i_OK + i_IK,
-        D(i_CK1) ~ rc0c1 - rc1c2,
-        D(i_CK2) ~ rc1c2 - rc2o,
-        D(i_OK) ~ rc2o - roi,
-        D(i_IK) ~ roi,
-        yinf ~ expit(-(V + 78.65) / 6.33),
-        tauy ~ 1 / (0.11885 * exp((V + 75) / 28.37) + 0.56236 * exp(-(V + 75) / 14.19)),
-        IfNa ~ gf * fNa * i_y * (vm - E_Na),
-        IfK ~ gf * fK * i_y * (vm - E_K),
-        If ~ IfNa + IfK,
-        D(i_y) ~ (yinf - i_y)/tauy
-    ], t; name)
+            E_Na ~ nernst(na_o, na_i),
+            E_K ~ nernst(k_o, k_i),
+            IK1 ~ gK1 * vk1 / (0.1653 + exp(0.0319 / mV * vk1)),
+            sinf ~ expit((V + 31.97156) / -4.64291),
+            rinf ~ expit((V - 3.55716) / 14.61299),
+            slowinf ~ sinf,
+            taur ~ inv(45.16 * exp(0.03577 * (V + 50)) + 98.9 * exp(-0.1 * (V + 38))),
+            taus ~ (0.35 * exp(-(((V + 70) / 15)^2)) + 0.035) - 26.9ms,
+            tausslow ~ (3.7 * exp(-(((V + 70) / 30)^2)) + 0.035) + 37.4ms,
+            Ito ~ gt * i_r * (f_is * i_s + (1 - f_is) * i_sslow) * (vm - E_K),
+            D(i_r) ~ (rinf - i_r) / taur,
+            D(i_s) ~ (sinf - i_s) / taus,
+            D(i_sslow) ~ (slowinf - i_sslow) / tausslow,
+            IKs ~ GKs * i_nKs^2 * (vm - E_K) * fracIKuravail * 2,
+            D(i_nKs) ~ (alphan / (alphan + betan) - i_nKs) / nKstau,
+            E_Kr ~ nernst(0.98 * k_o + 0.02 * na_o, 0.98 * k_i + 0.02 * na_i),
+            IKr ~ GKr * i_OK * (vm - E_Kr),
+            1 ~ CK0 + i_CK1 + i_CK2 + i_OK + i_IK,
+            D(i_CK1) ~ rc0c1 - rc1c2,
+            D(i_CK2) ~ rc1c2 - rc2o,
+            D(i_OK) ~ rc2o - roi,
+            D(i_IK) ~ roi,
+            yinf ~ expit(-(V + 78.65) / 6.33),
+            tauy ~ 1 / (0.11885 * exp((V + 75) / 28.37) + 0.56236 * exp(-(V + 75) / 14.19)),
+            IfNa ~ gf * fNa * i_y * (vm - E_Na),
+            IfK ~ gf * fK * i_y * (vm - E_K),
+            If ~ IfNa + IfK,
+            D(i_y) ~ (yinf - i_y) / tauy
+        ], t; name)
 end
- # Ito: Where does this come from? Perhaps here: https://modeldb.science/262081
+# Ito: Where does this come from? Perhaps here: https://modeldb.science/262081
diff --git a/src/ina.jl b/src/ina.jl
index ab33697..6efee8c 100644
--- a/src/ina.jl
+++ b/src/ina.jl
@@ -19,7 +19,7 @@ function get_ina_sys(nai, nao, vm; name=:inasys)
         E_Na(t)
     end
 
-    V = vm * Volt // mV # Convert voltage to mV
+    V = vm * inv(mV) # Convert voltage quantity to mV
     NatauhHI = 0.4537ms * expit((V + 10.66) / 11.1)
     NatauhLOW = 3.49ms / (0.135 * exp((V + 80) / -6.8) + 3.56 * exp(0.079V) + 3.1e5 * exp(0.35V))
     NataujHI = 11.63ms * (1 + exp(-0.1 * (V + 32))) / exp(-2.535e-7V)
diff --git a/src/isoproterenol.jl b/src/isoproterenol.jl
index 5fa3431..7d3450f 100644
--- a/src/isoproterenol.jl
+++ b/src/isoproterenol.jl
@@ -1,15 +1,7 @@
-# Beta adrenergic system activated by isoproterenol
-using ModelingToolkit
-using ModelingToolkit: t_nounits as t, D_nounits as D
-
-function get_bar_sys(;
-    ATP=5000μM,
-    ISO=0μM,
-    name=:bar_sys,
-    simplify=false
-)
+"Beta adrenergic system activated by isoproterenol"
+function get_bar_sys(ATP=5000μM, ISO=0μM; name=:bar_sys, simplify=false)
     @parameters begin
-        b1ARtot = 0.00528μM
+        b1ARtot = 5.28e-3μM
         Gstot = 3.83μM
         PDEtot = 22.85e-3μM
         PKItot = 0.18μM
@@ -53,11 +45,11 @@ function get_bar_sys(;
         k_PP_PDE = 1.5Hz                         # rate constant for PDE dephosphorylation by phosphatases
         kf_RC_cAMP = 1 / (μM * ms)               # Kd for PKA RC binding to cAMP
         kf_RCcAMP_cAMP = 1 / (μM * ms)           # Kd for PKA RC:cAMP binding to cAMP
-        kf_RcAMPcAMP_C = 4.375 / (μM * ms)       # Kd for PKA R:cAMPcAMP binding to C
+        kf_RcAMPcAMP_C = 4.375 / ms       # Kd for PKA R:cAMPcAMP binding to C
         kf_PKA_PKI = 1 / (μM * ms)               # Ki for PKA inhibition by PKI
         kr_RC_cAMP = 1.64 / ms                   # rate constant for PKA RC unbinding to cAMP
         kr_RCcAMP_cAMP = 9.14 / ms               # rate constant for PKA RC:cAMP unbinding to cAMP
-        kr_RcAMPcAMP_C = 1 / ms                  # rate constant for PKA R:cAMPcAMP unbinding to C
+        kr_RcAMPcAMP_C = 1 / (μM * ms)           # rate constant for PKA R:cAMPcAMP binding to C
         kr_PKA_PKI = 0.2Hz                       # reverse rate for PKA inhibition by PKI
         epsilon = 10                             # AKAP-mediated scaling factor
         k_PKA_I1 = 60Hz                          # rate constant for I-1 phosphorylation by type 1 PKA
diff --git a/src/istim.jl b/src/istim.jl
index 8e6b537..6b7bb1a 100644
--- a/src/istim.jl
+++ b/src/istim.jl
@@ -1,13 +1,11 @@
-using DiffEqCallbacks
-
-function build_stim_callbacks(starttime, endtime, period, duty, strength, sym; baseline=0, proposeddt=0.01)
+function build_stim_callbacks(sym, endtime; period = 1second, duty = 5e-4second,starttime=zero(endtime), strength=-80μAμF, baseline=0μAμF, proposeddt=0.01second)
     rise! = (integrator) -> begin
-        integrator[sym] = strength
+        integrator.ps[sym] = strength
         set_proposed_dt!(integrator, proposeddt)
     end
 
     fall! = (integrator) -> begin
-        integrator[sym] = baseline
+        integrator.ps[sym] = baseline
         set_proposed_dt!(integrator, proposeddt)
     end
 
diff --git a/src/utils.jl b/src/utils.jl
index a136460..c4d9f07 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1,41 +1,42 @@
 # Units
 const second = 1           # second
 const minute = 60second    # minute
-const ms = second//1000    # millisecond
-const Hz = 1 // second     # Herz
+const ms =  0.001second    # millisecond
+const Hz =  inv(second)    # Herz
 const kHz = 1000Hz         # kilohertz
 const metre = 1            # meter
-const cm = metre//100      # centimeter
+const cm = 0.01metre       # centimeter
 const cm² = cm^2           # square centimeter
-const μm = metre//10^6     # Micrometer
+const μm = 1E-6metre       # Micrometer
 const mL = cm^3            # milliliter = cubic centimeter
 const Liter = 1000mL       # liter
-const μL = Liter//10^6
-const pL = Liter//10^12
-const mM = 1
-const Molar = 1000mM       # molar (1000 since the SI units is mM)
-const μM = mM//10^3        # micromolar
-const nM = mM//10^6        # nanomolar
+const μL = 1E-6Liter
+const pL = 1E-12Liter
+const mM = 1               # the SI units is mM
+const Molar = 1000mM       # molar
+const μM = 0.001mM         # micromolar
+const nM = 1E-6mM          # nanomolar
 const Amp = 1              # ampere
-const mA = Amp//10^3         # milliampere
-const μA = Amp//10^6         # micropampere
-const Volt = 1               # volt
-const mV = Volt//10^3        # millivolt
-const mS = mA // Volt        # milliseimens
-const Farad = Amp * second // Volt
-const μF = Farad//10^6
-const T₀ = 310                  # Default temp (37C)
-const Faraday = 96485           # Faraday constant (columb / mol)
-const RGAS = 8314//1000         # Ideal gas constant (J/K⋅mol)
-const VT = RGAS * T₀ // Faraday # Thermal voltage (@37C), about 26.7 mV
-const iVT = inv(VT)             # Reciprocal of thermal voltage
-const μAμF = μA // μF           # Common unit for current density, normalized by capacitance
-const mSμF = ms // μF           # Common unit for conductance, normalized by capacitance
+const mA = 0.001Amp        # milliampere
+const μA = 1E-6Amp         # micropampere
+const Volt = 1             # volt
+const mV = 0.001Volt       # millivolt
+const mS = mA / Volt        # milliseimens
+const Farad = Amp * second / Volt
+const μF = 1E-6Farad
+const T₀ = 310                 # Default temp (37C)
+const Faraday = 96485          # Faraday constant (columb / mol)
+const RGAS = 8.314             # Ideal gas constant (J/K⋅mol)
+const VT = RGAS * T₀ / Faraday # Thermal voltage (@37C), about 26.7 mV
+const iVT = inv(VT)            # Reciprocal of thermal voltage
+const μAμF = μA / μF           # Common unit for current density, normalized by capacitance
+const mSμF = ms / μF           # Common unit for conductance, normalized by capacitance
 
 """
 Regular Hill/MM function
 """
 hil(x, k=one(x)) = x / (x + k)
+hil(x, k, n::Int) = hil(x^n, k^n)
 hil(x, k, n) = hil(NaNMath.pow(x, n), NaNMath.pow(k, n))
 
 """
@@ -56,7 +57,7 @@ exprel(x, em1=expm1(x)) = x / em1
 
 """Nernst potential"""
 nernst(x_out, x_in) = VT * NaNMath.log(x_out / x_in)
-nernst(x_out, x_in, z) = nernst(x_out, x_in) // z
+nernst(x_out, x_in, z) = nernst(x_out, x_in) / z
 
 "GHK flux equation"
 ghk(px, x_i, x_o, zvfrt, ezvfrtm1=expm1(zvfrt), z=1) = px * z * Faraday * ((ezvfrtm1 + 1) * x_i - x_o) * exprel(zvfrt, ezvfrtm1)
diff --git a/test.jl b/test.jl
index 9b9f01f..a37344a 100644
--- a/test.jl
+++ b/test.jl
@@ -1,20 +1,24 @@
-using CaMKIIModel
+using ProgressLogging
 using ModelingToolkit
-using OrdinaryDiffEq
-using BenchmarkTools
+using DifferentialEquations
 using Plots
+using CaMKIIModel
+Plots.default(lw=2)
 
 sys = build_neonatal_ecc_sys(simplify=true)
+tend = 1000.0
+prob = ODEProblem(sys, [], tend)
 
-prob = ODEProblem(sys, [], float(10))
-
-# TODO: fix unstable (ICaL?)
-sol = solve(prob, Rodas5P())
+@unpack Istim = sys
+callback = build_stim_callbacks(Istim, tend)
+alg = FBDF()
+@time sol = solve(prob, alg; callback, progress=true, abstol=1e-6, reltol=1e-6, maxiters=Int(1e8))
 
-sol[sys.ICaL]
 
-plot(sol, idxs=sys.PO1RyR)
+plot(sol, idxs=sys.vm, tspan=(90, 100))
+plot(sol, idxs=[sys.CaJSR, sys.CaNSR], tspan=(90, 100))
+plot(sol, idxs=[sys.Cai_sub_SR, sys.Cai_sub_SL, sys.Cai_mean], tspan=(99, 100))
+plot(sol, idxs=sys.Jtr)
 
-for (k, v) in zip(unknowns(sys), sol[end])
-    println(k, " => ", v)
-end
+u0 = sol[end]
+@time sol = solve(remake(prob, u0=u0), alg; progress=true, abstol=1e-9, reltol=1e-9)