Caffeine Effects - NRVM CaMKII model

Caffeine increase RyR opening sensitivity to luminal and subspace calcium. In this model, we decrease the mid saturation sub-SR calcium concentration for the opening rate.

using ModelingToolkit
using OrdinaryDiffEq, SteadyStateDiffEq, DiffEqCallbacks
using Plots
using CSV
using DataFrames
import Dates
using Model
using Model: second
Plots.default(lw=1.5)

@time "Build system" sys = Model.DEFAULT_SYS
tend = 500second
@time "Build problem" prob = ODEProblem(sys, [], tend)
stimstart = 100second
stimend = 300second
alg = KenCarp47()
function add_coffee_affect!(integrator)
    integrator.ps[sys.RyRsensitivity] = 10
end

Build system: 0.000006 seconds
Build problem: 52.361792 seconds (81.23 M allocations: 5.496 GiB, 3.59% gc time, 99.62% compilation time: 17% of which was recompilation)

add_coffee_affect! (generic function with 1 method)

@unpack Istim = sys
callback = build_stim_callbacks(Istim, stimend; period=1second, starttime=stimstart)

SciMLBase.CallbackSet{Tuple{}, Tuple{SciMLBase.DiscreteCallback{DiffEqCallbacks.PresetTimeFunction{StepRange{Int64, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#15#17"{Float64, Float64, Symbolics.Num}}, Model.var"#15#17"{Float64, Float64, Symbolics.Num}, DiffEqCallbacks.PresetTimeFunction{StepRange{Int64, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#15#17"{Float64, Float64, Symbolics.Num}}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}, SciMLBase.DiscreteCallback{DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#16#18"{Float64, Float64, Symbolics.Num}}, Model.var"#16#18"{Float64, Float64, Symbolics.Num}, DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#16#18"{Float64, Float64, Symbolics.Num}}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}}}((), (SciMLBase.DiscreteCallback{DiffEqCallbacks.PresetTimeFunction{StepRange{Int64, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#15#17"{Float64, Float64, Symbolics.Num}}, Model.var"#15#17"{Float64, Float64, Symbolics.Num}, DiffEqCallbacks.PresetTimeFunction{StepRange{Int64, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#15#17"{Float64, Float64, Symbolics.Num}}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}(DiffEqCallbacks.PresetTimeFunction{StepRange{Int64, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#15#17"{Float64, Float64, Symbolics.Num}}(100000:1000:300000, true, SciMLBase.INITIALIZE_DEFAULT, Model.var"#15#17"{Float64, Float64, Symbolics.Num}(-80.0, 0.5, Istim)), Model.var"#15#17"{Float64, Float64, Symbolics.Num}(-80.0, 0.5, Istim), DiffEqCallbacks.PresetTimeFunction{StepRange{Int64, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#15#17"{Float64, Float64, Symbolics.Num}}(100000:1000:300000, true, SciMLBase.INITIALIZE_DEFAULT, Model.var"#15#17"{Float64, Float64, Symbolics.Num}(-80.0, 0.5, Istim)), SciMLBase.FINALIZE_DEFAULT, Bool[1, 1], nothing, (), true), SciMLBase.DiscreteCallback{DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#16#18"{Float64, Float64, Symbolics.Num}}, Model.var"#16#18"{Float64, Float64, Symbolics.Num}, DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#16#18"{Float64, Float64, Symbolics.Num}}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}(DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#16#18"{Float64, Float64, Symbolics.Num}}(100000.5:1000.0:300000.5, true, SciMLBase.INITIALIZE_DEFAULT, Model.var"#16#18"{Float64, Float64, Symbolics.Num}(0.0, 0.5, Istim)), Model.var"#16#18"{Float64, Float64, Symbolics.Num}(0.0, 0.5, Istim), DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#16#18"{Float64, Float64, Symbolics.Num}}(100000.5:1000.0:300000.5, true, SciMLBase.INITIALIZE_DEFAULT, Model.var"#16#18"{Float64, Float64, Symbolics.Num}(0.0, 0.5, Istim)), SciMLBase.FINALIZE_DEFAULT, Bool[1, 1], nothing, (), true)))

Add caffeine at t = 200 econd

callback_caf = CallbackSet(build_stim_callbacks(Istim, stimend; period=1second, starttime=stimstart), PresetTimeCallback(200.0second, add_coffee_affect!))

SciMLBase.CallbackSet{Tuple{}, Tuple{SciMLBase.DiscreteCallback{DiffEqCallbacks.PresetTimeFunction{StepRange{Int64, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#15#17"{Float64, Float64, Symbolics.Num}}, Model.var"#15#17"{Float64, Float64, Symbolics.Num}, DiffEqCallbacks.PresetTimeFunction{StepRange{Int64, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#15#17"{Float64, Float64, Symbolics.Num}}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}, SciMLBase.DiscreteCallback{DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#16#18"{Float64, Float64, Symbolics.Num}}, Model.var"#16#18"{Float64, Float64, Symbolics.Num}, DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#16#18"{Float64, Float64, Symbolics.Num}}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}, SciMLBase.DiscreteCallback{DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(Main.var"##225".add_coffee_affect!)}, typeof(Main.var"##225".add_coffee_affect!), DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(Main.var"##225".add_coffee_affect!)}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}}}((), (SciMLBase.DiscreteCallback{DiffEqCallbacks.PresetTimeFunction{StepRange{Int64, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#15#17"{Float64, Float64, Symbolics.Num}}, Model.var"#15#17"{Float64, Float64, Symbolics.Num}, DiffEqCallbacks.PresetTimeFunction{StepRange{Int64, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#15#17"{Float64, Float64, Symbolics.Num}}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}(DiffEqCallbacks.PresetTimeFunction{StepRange{Int64, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#15#17"{Float64, Float64, Symbolics.Num}}(100000:1000:300000, true, SciMLBase.INITIALIZE_DEFAULT, Model.var"#15#17"{Float64, Float64, Symbolics.Num}(-80.0, 0.5, Istim)), Model.var"#15#17"{Float64, Float64, Symbolics.Num}(-80.0, 0.5, Istim), DiffEqCallbacks.PresetTimeFunction{StepRange{Int64, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#15#17"{Float64, Float64, Symbolics.Num}}(100000:1000:300000, true, SciMLBase.INITIALIZE_DEFAULT, Model.var"#15#17"{Float64, Float64, Symbolics.Num}(-80.0, 0.5, Istim)), SciMLBase.FINALIZE_DEFAULT, Bool[1, 1], nothing, (), true), SciMLBase.DiscreteCallback{DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#16#18"{Float64, Float64, Symbolics.Num}}, Model.var"#16#18"{Float64, Float64, Symbolics.Num}, DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#16#18"{Float64, Float64, Symbolics.Num}}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}(DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#16#18"{Float64, Float64, Symbolics.Num}}(100000.5:1000.0:300000.5, true, SciMLBase.INITIALIZE_DEFAULT, Model.var"#16#18"{Float64, Float64, Symbolics.Num}(0.0, 0.5, Istim)), Model.var"#16#18"{Float64, Float64, Symbolics.Num}(0.0, 0.5, Istim), DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Model.var"#16#18"{Float64, Float64, Symbolics.Num}}(100000.5:1000.0:300000.5, true, SciMLBase.INITIALIZE_DEFAULT, Model.var"#16#18"{Float64, Float64, Symbolics.Num}(0.0, 0.5, Istim)), SciMLBase.FINALIZE_DEFAULT, Bool[1, 1], nothing, (), true), SciMLBase.DiscreteCallback{DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(Main.var"##225".add_coffee_affect!)}, typeof(Main.var"##225".add_coffee_affect!), DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(Main.var"##225".add_coffee_affect!)}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}(DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(Main.var"##225".add_coffee_affect!)}([200000.0], true, SciMLBase.INITIALIZE_DEFAULT, Main.var"##225".add_coffee_affect!), Main.var"##225".add_coffee_affect!, DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(Main.var"##225".add_coffee_affect!)}([200000.0], true, SciMLBase.INITIALIZE_DEFAULT, Main.var"##225".add_coffee_affect!), SciMLBase.FINALIZE_DEFAULT, Bool[1, 1], nothing, (), true)))

Single-dose caffeine¶

@time sol = solve(prob, alg; callback)
@time sol_caf = solve(prob, alg; callback=callback_caf)

i = (sys.t / 1000, sys.vm)
plot(sol, idxs=i, title="Action potential", lab="Ctl", tspan=(198second, 205second))
plot!(sol_caf, idxs=i, lab="Caf", tspan=(198second, 205second), ylabel="Voltage (mV)", xlabel="Time (s)")

  6.128804 seconds (8.34 M allocations: 972.445 MiB, 2.88% gc time, 64.71% compilation time)
  4.806327 seconds (4.29 M allocations: 707.032 MiB, 1.65% gc time, 53.96% compilation time)

i = (sys.t / 1000, sys.Cai_sub_SR * 1000)
plot(sol, idxs=i, title="Calcium transient (During caffeine addition)", lab="Ctl", tspan=(198second, 205second))
plot!(sol_caf, idxs=i, tspan=(198second, 205second), lab="Caf", ylabel="Subspace calcium (nM)", xlabel="Time (s)")

i = (sys.t / 1000, sys.PO1RyR)
plot(sol, idxs=i, title="RyR open (During caffeine addition)", lab="Ctl", tspan=(198second, 205second))
plot!(sol_caf, idxs=i, tspan=(198second, 205second), lab="Caf", ylabel="Open probability", ylims=(0, 1), xlabel="Time (s)")

i = (sys.t / 1000, sys.Cai_sub_SR * 1000)
plot(sol, idxs=i, title="Calcium transient (After caffeine addition)", lab="Ctl", ylabel="Subspace calcium (nM)", tspan=(198second, 205second))
plot!(sol_caf, idxs=i, lab="Caf", xlabel="Time (s)", tspan=(198second, 205second))

i = (sys.t / 1000, sys.CaJSR)
plot(sol, idxs=i, title="SR Calcium (During caffeine addition)", lab="Ctl", ylabel="SR calcium (μM)", tspan=(198second, 205second))
plot!(sol_caf, idxs=i, tspan=(198second, 205second), lab="Caf", ylims=(0, 850), xlabel="Time (s)")

i = (sys.t / 1000, sys.Jrel)
plot(sol, idxs=sys.Jrel, title="Ca flux", lab="Ctl  (Jrel)", tspan=(198second, 205second))
plot!(sol_caf, idxs=sys.Jrel, lab="Caf (Jrel)", tspan=(198second, 205second), ylabel="μM/ms", xlabel="Time (s)")

i = (sys.t / 1000, sys.CaMKAct * 100)
plot(sol, idxs=i, title="Active CaMKII", lab="Ctl")
plot!(sol_caf, idxs=i, lab="Caf", ylabel="CaMKII activity (%)", xlabel="Time (s)")

Caffeine and electrophysiology¶

Add caffeine in the beginning of the simulation.
Add caffeine and nifedipine in the beginning of the simulation (nifedipine blocks 90% of L-type calcium channel).

@time "Build system" sys = Model.DEFAULT_SYS
tend = 205second
stimstart = 30second
stimend = 120second
alg = KenCarp47()
@unpack Istim = sys
callback = build_stim_callbacks(Istim, stimend; period=1second, starttime=stimstart)
@time prob = ODEProblem(sys, [], tend)
@time prob_caf = ODEProblem(sys, [sys.RyRsensitivity => 10], tend)
gCaL = prob.ps[sys.GCaL]
@time prob_nif_caf = ODEProblem(sys, [sys.RyRsensitivity => 10, sys.GCaL => 0.1 * gCaL], tend)

ssalg = DynamicSS(alg)
sprob_caf = SteadyStateProblem(prob_caf)
sssol = solve(sprob_caf, ssalg; abstol=1e-10, reltol=1e-10)
@time sol = solve(prob, alg; callback)
@time sol_caf = solve(remake(prob_caf, u0=sssol.u), alg; callback)
@time sol_nif_caf = solve(remake(prob_nif_caf, u0=sssol.u), alg; callback)

i = (sys.t / 1000, sys.vm)
tspan = (100second, 102second)
plot(sol, idxs=i, title="Action potential", lab="Ctl"; tspan)
plot!(sol_caf, idxs=i, lab="Caf"; tspan)
plot!(sol_nif_caf, idxs=i, lab="Caf + Nif", tspan=tspan, ylabel="Voltage (mV)", xlabel="Time (s)")

Build system: 0.000001 seconds
  0.201712 seconds (856.75 k allocations: 43.577 MiB, 8.13% gc time, 27.41% compilation time)
  0.296233 seconds (892.69 k allocations: 46.762 MiB, 58.32% compilation time)
  0.510226 seconds (1.11 M allocations: 63.006 MiB, 3.75% gc time, 76.69% compilation time)
  1.041410 seconds (507.42 k allocations: 223.430 MiB, 4.10% gc time)
  1.119317 seconds (696.49 k allocations: 236.423 MiB, 1.12% gc time, 10.98% compilation time)
  0.582159 seconds (310.68 k allocations: 129.447 MiB)

savefig("caf-ap.pdf")

"/home/runner/work/camkii-cardiomyocyte-model/camkii-cardiomyocyte-model/docs/caf-ap.pdf"

i = (sys.t / 1000, sys.Cai_sub_SR * 1000)
tspan = (100second, 102second)
plot(sol, idxs=i, title="Calcium transient", lab="Ctl"; tspan)
plot!(sol_caf, idxs=i, lab="Caf"; tspan)
plot!(sol_nif_caf, idxs=i, lab="Caf + Nif", ylabel="Subspace calcium (nM)", xlabel="Time (s)"; tspan)

savefig("caf-cat.pdf")

"/home/runner/work/camkii-cardiomyocyte-model/camkii-cardiomyocyte-model/docs/caf-cat.pdf"

i = (sys.t / 1000, sys.CaJSR)
tspan = (100second, 102second)
plot(sol, idxs=i, title="SR Calcium", lab="Ctl", ylabel="SR calcium (μM)"; tspan)
plot!(sol_caf, idxs=i, lab="Caf"; tspan)
plot!(sol_nif_caf, idxs=i, lab="Caf + Nif", ylims=(0, 850), xlabel="Time (s)"; tspan)

CaMKII activities¶

Simulations¶

i = (sys.t / 1000, sys.CaMKAct * 100)
plot(sol, idxs=i, title="Simulation", lab="Ctl")
plot!(sol_caf, idxs=i, lab="Caf")
plot!(sol_nif_caf, idxs=i, lab="Caf + Nif", ylabel="CaMKII active fraction (%)", xlabel="Time (s)")

savefig("caf-camkact.pdf")

"/home/runner/work/camkii-cardiomyocyte-model/camkii-cardiomyocyte-model/docs/caf-camkact.pdf"

Experiments¶

chemicaldf = CSV.read(joinpath(@__DIR__, "data/CaMKAR-chemical.csv"), DataFrame)
ts = Dates.value.(chemicaldf[!, "Time"]) ./ 10^9
ctl = chemicaldf[!, "Ctrl Mean"]
ctl_error = chemicaldf[!, "Ctrl SD"] ./ sqrt.(chemicaldf[!, "Ctrl N"])

caf = chemicaldf[!, "caffeine 20mM Mean"]
caf_error = chemicaldf[!, "caffeine 20mM SD"] ./ sqrt.(chemicaldf[!, "caffeine 20mM N"])

plot(ts, ctl, yerr=ctl_error, lab="Control", color=:blue, markerstrokecolor=:blue)
plot!(ts, caf, yerr=caf_error, lab="Caffeine 20mM", color=:red, markerstrokecolor=:red)
plot!(xlabel="Time (s)", ylabel="CaMKII activity (A.U.)", title= "Experiment")

savefig("caf-exp.pdf")

"/home/runner/work/camkii-cardiomyocyte-model/camkii-cardiomyocyte-model/docs/caf-exp.pdf"

This notebook was generated using Literate.jl.