diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 85b1656..806a74a 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -91,7 +91,10 @@ jobs: key: ${{ steps.cache-julia.outputs.cache-primary-key }} - name: Run notebooks if: ${{ steps.cache-nb.outputs.cache-hit != 'true' }} - run: julia --project=@. --color=yes -p ${{ env.LITERATE_PROC }} ci.jl + run: | + mkdir -p ${{ env.NBCACHE }}/docs + cp - r docs/data ${{ env.NBCACHE }}/docs/ + julia --project=@. --color=yes -p ${{ env.LITERATE_PROC }} ci.jl - name: Copy back built notebooks run: cp --verbose -rf ${{ env.NBCACHE }}/docs/* docs/ - name: Build website diff --git a/docs/bcl.ipynb b/docs/bcl.ipynb deleted file mode 100644 index 1550d2e..0000000 --- a/docs/bcl.ipynb +++ /dev/null @@ -1,392 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Pacing frequency" - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [], - "source": [ - "using ModelingToolkit\n", - "using OrdinaryDiffEq\n", - "using Plots\n", - "using CSV\n", - "using DataFrames\n", - "using CaMKIIModel\n", - "Plots.default(lw=1.5, fmt=:png)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Setup the ODE system\n", - "\n", - "Electrical stimulation starts at `t`=100 seconds and ends at `t`=300 seconds." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "sys = build_neonatal_ecc_sys(simplify=true, reduce_iso=true)\n", - "tend = 500.0\n", - "prob = ODEProblem(sys, [], tend)\n", - "stimstart = 100.0\n", - "stimend = 300.0\n", - "@unpack Istim = sys\n", - "alg = FBDF()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 1Hz" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "callback = build_stim_callbacks(Istim, stimend; period=1, starttime=stimstart)\n", - "sol = solve(prob, alg; callback, abstol=1e-6, reltol=1e-6)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot(sol, idxs=sys.vm*1000, title=\"Action potential\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot(sol, idxs=sys.vm*1000, title=\"Action potential\", tspan=(299, 300))" - ] - }, - { - "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 transient\")" - ] - }, - { - "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, stimend; period=1/2, starttime=stimstart)\n", - "sol2 = solve(prob, alg; callback, abstol=1e-6, reltol=1e-6)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot(sol2, idxs=sys.vm*1000, title=\"Action potential\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot(sol2, idxs=sys.vm*1000, title=\"Action potential\", tspan=(299, 300))" - ] - }, - { - "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 transient\")" - ] - }, - { - "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, stimend; period=1/3, starttime=stimstart)\n", - "sol3 = solve(prob, alg; callback, abstol=1e-6, reltol=1e-6)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot(sol3, idxs=sys.vm*1000, title=\"Action potential\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot(sol3, idxs=sys.vm*1000, title=\"Action potential\", tspan=(299, 300))" - ] - }, - { - "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 transient\")" - ] - }, - { - "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\n", - "Action potential" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot(sol, idxs=sys.vm*1000, title=\"Action potential\", lab=\"1Hz\")\n", - "plot!(sol2, idxs=sys.vm*1000, lab=\"2Hz\")\n", - "plot!(sol3, idxs=sys.vm*1000, lab=\"3Hz\", tspan=(299, 300), xlabel=\"Time (sec.)\", ylabel=\"Voltage (mV)\")" - ] - }, - { - "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, 1.0))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Data fitting\n", - "### Pacing duration and CaMKII activity" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "durationdf = CSV.read(\"data/CaMKAR-duration.csv\", DataFrame)\n", - "ts = durationdf[!, \"Time(sec)\"]\n", - "fifteen = durationdf[!, \"1Hz 15sec (Mean)\"]\n", - "fifteen_error = durationdf[!, \"1Hz 15sec (SD)\"] ./ sqrt.(durationdf[!, \"1Hz 15sec (N)\"])\n", - "thirty = durationdf[!, \"1Hz 30sec (Mean)\"]\n", - "thirty_error = durationdf[!, \"1Hz 30sec (SD)\"] ./ sqrt.(durationdf[!, \"1Hz 30sec (N)\"])\n", - "sixty = durationdf[!, \"1Hz 60sec (Mean)\"]\n", - "sixty_error = durationdf[!, \"1Hz 60sec (SD)\"] ./ sqrt.(durationdf[!, \"1Hz 60sec (N)\"])\n", - "ninety = durationdf[!, \"1Hz 90sec (Mean)\"]\n", - "ninety_error = durationdf[!, \"1Hz 90sec (SD)\"] ./ sqrt.(durationdf[!, \"1Hz 90sec (N)\"])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot(ts, fifteen, yerr=fifteen_error, lab=\"15 sec\", color = :blue, markerstrokecolor=:blue)\n", - "plot!(ts, thirty, yerr=thirty_error, lab=\"30 sec\", color = :red, markerstrokecolor=:red)\n", - "plot!(ts, sixty, yerr=sixty_error, lab=\"60 sec\", color = :orange, markerstrokecolor=:orange)\n", - "plot!(ts, ninety, yerr=ninety_error, lab=\"90 sec\", color = :green, markerstrokecolor=:green)\n", - "plot!(title=\"Pacing duration\", xlabel=\"Time (sec.)\", ylabel=\"CaMKII activity (AU)\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Resolution is reduced to 1Hz." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "stimstart=30.0\n", - "callback15 = build_stim_callbacks(Istim, stimstart+15; period=1,starttime=stimstart)\n", - "sol15 = solve(prob, alg; callback=callback15, abstol=1e-6, reltol=1e-6)\n", - "callback30 = build_stim_callbacks(Istim, stimstart+30; period=1, starttime=stimstart)\n", - "sol30 = solve(prob, alg; callback=callback30, abstol=1e-6, reltol=1e-6)\n", - "callback60 = build_stim_callbacks(Istim, stimstart+60; period=1, starttime=stimstart)\n", - "sol60 = solve(prob, alg; callback=callback60, abstol=1e-6, reltol=1e-6)\n", - "callback90 = build_stim_callbacks(Istim, stimstart+90; period=1, starttime=stimstart)\n", - "sol90 = solve(prob, alg; callback=callback90, abstol=1e-6, reltol=1e-6)\n", - "idx = sys.CaMKAct" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot(sol15, idxs=idx, tspan=(0, 205), lab=\"15 sec\", color = :blue)\n", - "plot!(sol30, idxs=idx, tspan=(0, 205), lab=\"30 sec\", color = :red)\n", - "plot!(sol60, idxs=idx, tspan=(0, 205), lab=\"60 sec\", color = :orange)\n", - "plot!(sol90, idxs=idx, tspan=(0, 205), lab=\"90 sec\", color = :green)\n", - "plot!(title=\"Pacing duration\", xlabel=\"Time (sec.)\", ylabel=\"CaMKII activity (AU)\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Pacing frequency and CaMKII activity" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "freqdf = CSV.read(\"data/CaMKAR-freq.csv\", DataFrame)\n", - "ts = 0:5:205\n", - "onehz = freqdf[!, \"1Hz (Mean)\"]\n", - "onehz_error = freqdf[!, \"1Hz (SD)\"] ./ sqrt.(freqdf[!, \"1Hz (N)\"])\n", - "twohz = freqdf[!, \"2Hz (Mean)\"]\n", - "twohz_error = freqdf[!, \"2Hz (SD)\"] ./ sqrt.(freqdf[!, \"2Hz (N)\"])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot(ts, onehz, yerr=onehz_error, lab=\"1 Hz\", color = :blue, markerstrokecolor=:blue)\n", - "plot!(ts, twohz, yerr=twohz_error, lab=\"2 Hz\", color = :red, markerstrokecolor=:red)\n", - "plot!(title=\"Pacing frequency\", xlabel=\"Time (sec.)\", ylabel=\"CaMKII activity (AU)\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tend = 205.0\n", - "prob = ODEProblem(sys, [], tend)\n", - "stimstart = 30.0\n", - "stimend = 120.0\n", - "callback = build_stim_callbacks(Istim, stimend; period=1, starttime=stimstart)\n", - "sol1 = solve(prob, alg; callback, abstol=1e-6, reltol=1e-6)\n", - "\n", - "callback2 = build_stim_callbacks(Istim, stimend; period=0.5, starttime=stimstart)\n", - "sol2 = solve(prob, alg; callback=callback2, abstol=1e-6, reltol=1e-6)\n", - "idx = sys.CaMKAct" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot(sol1, idxs=idx, lab=\"1 Hz\", color = :blue)\n", - "plot!(sol2, idxs=idx, lab=\"2 Hz\", color = :red)\n", - "plot!(title=\"Pacing frequency\", xlabel=\"Time (sec.)\", ylabel=\"CaMKII activity (AU)\", ylims=(0.0, 0.9))" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Julia 1.11.1", - "language": "julia", - "name": "julia-1.11" - }, - "language_info": { - "file_extension": ".jl", - "mimetype": "application/julia", - "name": "julia", - "version": "1.11.1" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/docs/bcl.jl b/docs/bcl.jl new file mode 100644 index 0000000..faab5c2 --- /dev/null +++ b/docs/bcl.jl @@ -0,0 +1,141 @@ +# # Pacing frequency +using ModelingToolkit +using OrdinaryDiffEq +using Plots +using CSV +using DataFrames +using CaMKIIModel +Plots.default(lw=1.5) + +# ## Setup the ODE system +# Electrical stimulation starts at `t`=100 seconds and ends at `t`=300 seconds. +sys = build_neonatal_ecc_sys(simplify=true, reduce_iso=true) +tend = 500.0 +prob = ODEProblem(sys, [], tend) +stimstart = 100.0 +stimend = 300.0 +@unpack Istim = sys +alg = FBDF() + +# ## 1Hz +callback = build_stim_callbacks(Istim, stimend; period=1, starttime=stimstart) +sol = solve(prob, alg; callback, abstol=1e-6, reltol=1e-6) + +#--- +plot(sol, idxs=sys.vm*1000, title="Action potential") + +#--- +plot(sol, idxs=sys.vm*1000, title="Action potential", tspan=(299, 300)) + +#--- +plot(sol, idxs=[sys.Cai_sub_SR, sys.Cai_sub_SL, sys.Cai_mean], tspan=(299, 300), title="Calcium transient") + +#--- +plot(sol, idxs=sys.CaMKAct, title="Active CaMKII") + +# ## 2Hz +callback = build_stim_callbacks(Istim, stimend; period=1/2, starttime=stimstart) +sol2 = solve(prob, alg; callback, abstol=1e-6, reltol=1e-6) + +#--- +plot(sol2, idxs=sys.vm*1000, title="Action potential") + +#--- +plot(sol2, idxs=sys.vm*1000, title="Action potential", tspan=(299, 300)) + +#--- +plot(sol2, idxs=[sys.Cai_sub_SR, sys.Cai_sub_SL, sys.Cai_mean], tspan=(299, 300), title="Calcium transient") + +#--- +plot(sol2, idxs=sys.CaMKAct, title="Active CaMKII") + +# ## 3Hz +callback = build_stim_callbacks(Istim, stimend; period=1/3, starttime=stimstart) +sol3 = solve(prob, alg; callback, abstol=1e-6, reltol=1e-6) + +#--- +plot(sol3, idxs=sys.vm*1000, title="Action potential") + +#--- +plot(sol3, idxs=sys.vm*1000, title="Action potential", tspan=(299, 300)) + +#--- +plot(sol3, idxs=[sys.Cai_sub_SR, sys.Cai_sub_SL, sys.Cai_mean], tspan=(299, 300), title="Calcium transient") + +#--- +plot(sol3, idxs=sys.CaMKAct, title="Active CaMKII") + +# ## Comparing 1-3 Hz +plot(sol, idxs=sys.vm*1000, title="Action potential", lab="1Hz") +plot!(sol2, idxs=sys.vm*1000, lab="2Hz") +plot!(sol3, idxs=sys.vm*1000, lab="3Hz", tspan=(299, 300), xlabel="Time (sec.)", ylabel="Voltage (mV)") + +#--- +plot(sol, idxs=sys.CaMKAct, title="CaMKII", lab="1Hz") +plot!(sol2, idxs=sys.CaMKAct, lab="2Hz") +plot!(sol3, idxs=sys.CaMKAct, lab="3Hz", ylim=(0.0, 1.0), xlabel="Time (sec.)", ylabel="Act. fraction (AU)") + +# ## Data fitting +### Pacing duration and CaMKII activity +durationdf = CSV.read(joinpath(@__DIR__, "data/CaMKAR-duration.csv"), DataFrame) +ts = durationdf[!, "Time(sec)"] +fifteen = durationdf[!, "1Hz 15sec (Mean)"] +fifteen_error = durationdf[!, "1Hz 15sec (SD)"] ./ sqrt.(durationdf[!, "1Hz 15sec (N)"]) +thirty = durationdf[!, "1Hz 30sec (Mean)"] +thirty_error = durationdf[!, "1Hz 30sec (SD)"] ./ sqrt.(durationdf[!, "1Hz 30sec (N)"]) +sixty = durationdf[!, "1Hz 60sec (Mean)"] +sixty_error = durationdf[!, "1Hz 60sec (SD)"] ./ sqrt.(durationdf[!, "1Hz 60sec (N)"]) +ninety = durationdf[!, "1Hz 90sec (Mean)"] +ninety_error = durationdf[!, "1Hz 90sec (SD)"] ./ sqrt.(durationdf[!, "1Hz 90sec (N)"]) + +plot(ts, fifteen, yerr=fifteen_error, lab="15 sec", color = :blue, markerstrokecolor=:blue) +plot!(ts, thirty, yerr=thirty_error, lab="30 sec", color = :red, markerstrokecolor=:red) +plot!(ts, sixty, yerr=sixty_error, lab="60 sec", color = :orange, markerstrokecolor=:orange) +plot!(ts, ninety, yerr=ninety_error, lab="90 sec", color = :green, markerstrokecolor=:green) +plot!(title="Pacing duration", xlabel="Time (sec.)", ylabel="CaMKII activity (AU)") + +# Simulation resolution is reduced to 1Hz. +stimstart=30.0 +callback15 = build_stim_callbacks(Istim, stimstart+15; period=1,starttime=stimstart) +sol15 = solve(prob, alg; callback=callback15, abstol=1e-6, reltol=1e-6) +callback30 = build_stim_callbacks(Istim, stimstart+30; period=1, starttime=stimstart) +sol30 = solve(prob, alg; callback=callback30, abstol=1e-6, reltol=1e-6) +callback60 = build_stim_callbacks(Istim, stimstart+60; period=1, starttime=stimstart) +sol60 = solve(prob, alg; callback=callback60, abstol=1e-6, reltol=1e-6) +callback90 = build_stim_callbacks(Istim, stimstart+90; period=1, starttime=stimstart) +sol90 = solve(prob, alg; callback=callback90, abstol=1e-6, reltol=1e-6) +idx = sys.CaMKAct + +plot(sol15, idxs=idx, tspan=(0, 205), lab="15 sec", color = :blue) +plot!(sol30, idxs=idx, tspan=(0, 205), lab="30 sec", color = :red) +plot!(sol60, idxs=idx, tspan=(0, 205), lab="60 sec", color = :orange) +plot!(sol90, idxs=idx, tspan=(0, 205), lab="90 sec", color = :green) +plot!(title="Pacing duration", xlabel="Time (sec.)", ylabel="CaMKII activity (AU)") + +# ### Pacing frequency and CaMKII activity +freqdf = CSV.read(joinpath(@__DIR__,"data/CaMKAR-freq.csv"), DataFrame) +ts = 0:5:205 +onehz = freqdf[!, "1Hz (Mean)"] +onehz_error = freqdf[!, "1Hz (SD)"] ./ sqrt.(freqdf[!, "1Hz (N)"]) +twohz = freqdf[!, "2Hz (Mean)"] +twohz_error = freqdf[!, "2Hz (SD)"] ./ sqrt.(freqdf[!, "2Hz (N)"]) + +plot(ts, onehz, yerr=onehz_error, lab="1 Hz", color = :blue, markerstrokecolor=:blue) +plot!(ts, twohz, yerr=twohz_error, lab="2 Hz", color = :red, markerstrokecolor=:red) +plot!(title="Pacing frequency", xlabel="Time (sec.)", ylabel="CaMKII activity (AU)") + +#--- +tend = 205.0 +prob = ODEProblem(sys, [], tend) +stimstart = 30.0 +stimend = 120.0 +callback = build_stim_callbacks(Istim, stimend; period=1, starttime=stimstart) +sol1 = solve(prob, alg; callback, abstol=1e-6, reltol=1e-6) + +callback2 = build_stim_callbacks(Istim, stimend; period=0.5, starttime=stimstart) +sol2 = solve(prob, alg; callback=callback2, abstol=1e-6, reltol=1e-6) +idx = sys.CaMKAct + +plot(sol1, idxs=idx, lab="1 Hz", color = :blue) +plot!(sol2, idxs=idx, lab="2 Hz", color = :red) +plot!(title="Pacing frequency", xlabel="Time (sec.)", ylabel="CaMKII activity (AU)", ylims=(0.0, 0.9))