Pacing durations - NRVM CaMKII model

Experiments vs simulations.

using ModelingToolkit
using ModelingToolkit: t_nounits as t
using OrdinaryDiffEq, SteadyStateDiffEq, DiffEqCallbacks
using Plots
using CSV
using DataFrames
using CurveFit
using Model
using Model: second, Hz
Plots.default(lw=1.5)

Experiments¶

30 sec resting + N sec 1Hz pacing + resting.

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)"] .+ 0.25
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)"] .- 0.25
ninety_error = durationdf[!, "1Hz 90sec (SD)"] ./ sqrt.(durationdf[!, "1Hz 90sec (N)"])

42-element Vector{Float64}:
 0.07322025867393339
 0.06288123803690504
 0.06438696849626417
 0.07094952584857149
 0.0722214908336562
 0.07401047293769683
 0.06323601113333674
 0.07088528768397788
 0.06905459448562148
 0.06667878213772092
 ⋮
 0.06219021582828689
 0.06524406587899219
 0.05082442404251369
 0.06307382411198
 0.07385218047440929
 0.062460484843557096
 0.03877222159045237
 0.05170103726195504
 0.042024632430084405

30 sec timeseries +0.25 and 90 sec timeseries -0.25 for a consistent baseline before pacing.

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="Experiment", xlabel="Time (s)", ylabel="CaMKII activity (AU)")

savefig("pacing-duration-exp.pdf")

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

Simulation¶

@time "Build system" sys = Model.DEFAULT_SYS
tend = 210second
@time "Build problem" prob = ODEProblem(sys, [], (0second, tend))
alg = KenCarp47()

Build system: 0.000003 seconds
Build problem: 41.719827 seconds (50.77 M allocations: 3.518 GiB, 3.37% gc time, 99.53% compilation time: 20% of which was recompilation)

KenCarp47(; linsolve = nothing, nlsolve = OrdinaryDiffEqNonlinearSolve.NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Nothing}(1//100, 10, 1//5, 1//5, false, true, nothing), precs = DEFAULT_PRECS, smooth_est = true, extrapolant = linear, controller = PI, autodiff = ADTypes.AutoForwardDiff(),)

stimstart = 30second
@unpack Istim, CaMKAct = sys
idxs = (t / 1000, CaMKAct * 100)

sols = map((15, 30, 60, 90)) do dur
    cb = build_stim_callbacks(Istim, stimstart + dur * second; period=1second, starttime=stimstart)
    sol = solve(prob, alg; callback=cb)
end

plot(title="Simulation")

for (sol, dur, color) in zip(sols, (15, 30, 60, 90), (:blue, :red, :orange, :green))
    plot!(sol, idxs=idxs, tspan=(0second, 205second), lab="$dur sec", color=color)
end

plot!(xlabel="Time (s)", ylabel="CaMKII activity (%)")

savefig("pacing-duration-sim.pdf")

"/home/runner/work/camkii-cardiomyocyte-model/camkii-cardiomyocyte-model/docs/pacing-duration-sim.pdf"

Decay rates¶

Fit data from experiments and simulations against an exponential decay model.

Record 50 seconds after pacing ends.

ts = collect(range(0.0, stop=50.0, step=5.0)) ## in seconds

ydata_15 = fifteen[10:20]
ydata_30 = thirty[13:23]
ydata_60 = sixty[19:29]
ydata_90 = ninety[25:35]

11-element Vector{Float64}:
 13.64164293
 13.48430494
 13.22906779
 13.12909541
 13.04527631
 12.89725458
 12.88348484
 12.81471802
 12.83336721
 12.74340616
 12.79848118

Simulation points

ysim_15 = sols[1](stimstart+15second:5second:stimstart+15second+50second ; idxs=sys.CaMKAct * 100).u
ysim_30 = sols[2](stimstart+30second:5second:stimstart+30second+50second ; idxs=sys.CaMKAct * 100).u
ysim_60 = sols[3](stimstart+60second:5second:stimstart+60second+50second ; idxs=sys.CaMKAct * 100).u
ysim_90 = sols[4](stimstart+90second:5second:stimstart+90second+50second ; idxs=sys.CaMKAct * 100).u

11-element Vector{Float64}:
 45.37316999714231
 36.201947849181835
 28.221164693040823
 22.580537576734628
 18.5135337345181
 15.493970521820849
 13.184186182129684
 11.36883273114771
  9.908671996559523
  8.711218173726776
  7.713718272731036

Fit

fit_15 = solve(CurveFitProblem(ts, ydata_15), ExpSumFitAlgorithm(n=1, withconst=true))
fit_30 = solve(CurveFitProblem(ts, ydata_30), ExpSumFitAlgorithm(n=1, withconst=true))
fit_60 = solve(CurveFitProblem(ts, ydata_60), ExpSumFitAlgorithm(n=1, withconst=true))
fit_90 = solve(CurveFitProblem(ts, ydata_90), ExpSumFitAlgorithm(n=1, withconst=true))

fit_sim_15 = solve(CurveFitProblem(ts, ysim_15), ExpSumFitAlgorithm(n=1, withconst=true))
fit_sim_30 = solve(CurveFitProblem(ts, ysim_30), ExpSumFitAlgorithm(n=1, withconst=true))
fit_sim_60 = solve(CurveFitProblem(ts, ysim_60), ExpSumFitAlgorithm(n=1, withconst=true))
fit_sim_90 = solve(CurveFitProblem(ts, ysim_90), ExpSumFitAlgorithm(n=1, withconst=true))

retcode: Success
alg: ExpSumFitAlgorithm
residuals mean: -2.2608177956003187e-15
u: [5.396425283137034, 40.06764162289509, -0.055316349464354105]

Fitting results (experiments)¶

p1 = plot(ts, ydata_15, label="Exp 15 sec")
plot!(p1, ts, predict(fit_15), label="Fit", linestyle=:dash)
p2 = plot(ts, ydata_30, label="Exp 30 sec")
plot!(p2, ts, predict(fit_30), label="Fit", linestyle=:dash)
p3 = plot(ts, ydata_60, label="Exp 60 sec")
plot!(p3, ts, predict(fit_60), label="Fit", linestyle=:dash)
p4 = plot(ts, ydata_90, label="Exp 90 sec")
plot!(p4, ts, predict(fit_90), label="Fit", linestyle=:dash)
plot(p1, p2, p3, p4, layout=(2,2), xlabel="Time (s)", ylabel="CaMKII activity (AU)")

Fitting results (simulations)¶

p1s = plot(ts, ysim_15, label="Sim 15 sec")
plot!(p1s, ts, predict(fit_sim_15), label="Fit", linestyle=:dash)
p2s = plot(ts, ysim_30, label="Sim 30 sec")
plot!(p2s, ts, predict(fit_sim_30), label="Fit", linestyle=:dash)
p3s = plot(ts, ysim_60, label="Sim 60 sec")
plot!(p3s, ts, predict(fit_sim_60), label="Fit", linestyle=:dash)
p4s = plot(ts, ysim_90, label="Sim 90 sec")
plot!(p4s, ts, predict(fit_sim_90), label="Fit", linestyle=:dash)
plot(p1s, p2s, p3s, p4s, layout=(2,2), xlabel="Time (s)", ylabel="CaMKII activity (%)")

Decay time scales (tau)¶

tau_exp_15 = inv(-fit_15.u.λ[])
tau_exp_30 = inv(-fit_30.u.λ[])
tau_exp_60 = inv(-fit_60.u.λ[])
tau_exp_90 = inv(-fit_90.u.λ[])
tau_sim_15 = inv(-fit_sim_15.u.λ[])
tau_sim_30 = inv(-fit_sim_30.u.λ[])
tau_sim_60 = inv(-fit_sim_60.u.λ[])
tau_sim_90 = inv(-fit_sim_90.u.λ[])

println("The time scales for experiments: ")
for (tau, dur) in zip((tau_exp_15, tau_exp_30, tau_exp_60, tau_exp_90), (15, 30, 60, 90))
    println("$dur sec pacing is $(round(tau; digits=2)) seconds.")
end

println("The time scales for simulations: ")
for (tau, dur) in zip((tau_sim_15, tau_sim_30, tau_sim_60, tau_sim_90), (15, 30, 60, 90))
    println("$dur sec pacing is $(round(tau; digits=2)) seconds.")
end

The time scales for experiments: 
15 sec pacing is 16.48 seconds.
30 sec pacing is 16.73 seconds.
60 sec pacing is 17.65 seconds.
90 sec pacing is 18.08 seconds.
The time scales for simulations: 
15 sec pacing is 16.35 seconds.
30 sec pacing is 17.16 seconds.
60 sec pacing is 17.68 seconds.
90 sec pacing is 18.08 seconds.

Plot pacing time vs decay time scale

pacing_durations = [15.0, 30.0, 60.0, 90.0]
tau_experiments = [tau_exp_15, tau_exp_30, tau_exp_60, tau_exp_90]
tau_simulations = [tau_sim_15, tau_sim_30, tau_sim_60, tau_sim_90]
plot(pacing_durations, tau_experiments, label="Experiments", marker=:circle, color=:blue)
plot!(pacing_durations, tau_simulations, label="Simulations", marker=:square, color=:red)
plot!(title="Pacing Duration vs Decay Time Scale", xlabel="Pacing Duration (s)", ylabel="Decay Time Scale (s)")

Phosphorylated fraction¶

idxs = (sys.t / 1000, (sys.CaMKP + sys.CaMKA + sys.CaMKA2) * 100)
tspan = (0second, 205second)
plot(sols[1], idxs=idxs, tspan=tspan, lab="15 sec", color=:blue)
plot!(sols[2], idxs=idxs, tspan=tspan, lab="30 sec", color=:red)
plot!(sols[3], idxs=idxs, tspan=tspan, lab="60 sec", color=:orange)
plot!(sols[4], idxs=idxs, tspan=tspan, lab="90 sec", color=:green)
plot!(title="Simulation", xlabel="Time (s)", ylabel="Phosphorylated CaMKII (%)")

This notebook was generated using Literate.jl.