ROS effects - NRVM CaMKII model

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

Setup model¶

@time "Build system" sys = Model.DEFAULT_SYS
@unpack Istim = sys
tend = 205second
stimstart = 30second
stimend = 120second
@time "Build problem" prob = ODEProblem(sys, [], tend)
callback = build_stim_callbacks(Istim, stimend; period=1second, starttime=stimstart)
alg = KenCarp47()

Build system: 0.000004 seconds
Build problem: 41.256356 seconds (50.78 M allocations: 3.518 GiB, 2.97% gc time, 99.46% 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(),)

Comparisons¶

@time "Solve problem" sol = solve(prob, alg; callback)

# ROS (H2O2) 0.1uM
prob2 = remake(prob, p=[sys.ROS => 0.1μM])
@time "Solve problem" sol2 = solve(prob2, alg; callback)

# ROS (H2O2) 0.5uM
prob3 = remake(prob, p=[sys.ROS => 0.5μM])
@time "Solve problem" sol3 = solve(prob3, alg; callback);

Solve problem: 4.791011 seconds (7.77 M allocations: 723.000 MiB, 2.90% gc time, 78.41% compilation time)
Solve problem: 1.011221 seconds (502.70 k allocations: 219.661 MiB, 1.69% gc time)
Solve problem: 1.028016 seconds (501.85 k allocations: 218.946 MiB, 1.11% gc time)

Comparisons¶

i = (sys.t / 1000, sys.CaMKAct * 100)
plot(sol, idxs=i, lab="ROS (-)")
plot!(sol2, idxs=i, lab="ROS 0.1uM")
plot!(sol3, idxs=i, lab="ROS 0.5uM")
plot!(xlabel="Time (s)", ylabel="Active fraction (%)", title="Simulation")

savefig("ros-camkii.pdf")

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

Oxidized fraction

i = (sys.t / 1000, 100 * (sys.CaMKBOX + sys.CaMKPOX + sys.CaMKAOX + sys.CaMKOX ))
plot(sol, idxs=i, lab="ROS (-)")
plot!(sol2, idxs=i, lab="ROS 0.1uM")
plot!(sol3, idxs=i, lab="ROS 0.5uM")
plot!(xlabel="Time (s)", ylabel="Oxidized fraction (%)", title="Simulation")

savefig("ros-camkiiox.pdf")

"/home/runner/work/camkii-cardiomyocyte-model/camkii-cardiomyocyte-model/docs/ros-camkiiox.pdf"

Autophosphorylated fraction

i = (sys.t / 1000, 100 * (sys.CaMKP + sys.CaMKA + sys.CaMKA2))
plot(sol, idxs=i, lab="ROS (-)")
plot!(sol2, idxs=i, lab="ROS 0.1uM")
plot!(sol3, idxs=i, lab="ROS 0.5uM")
plot!(xlabel="Time (s)", ylabel="Phosphorylated fraction (%)", title="Simulation")

savefig("ros-camkiip.pdf")

"/home/runner/work/camkii-cardiomyocyte-model/camkii-cardiomyocyte-model/docs/ros-camkiip.pdf"

Experimental data¶

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"])

42-element Vector{Float64}:
 0.02305850621757014
 0.023225557529411765
 0.023536003837902393
 0.022995786588235295
 0.023056161133316953
 0.02236794743822152
 0.02448575734900888
 0.02631837493460871
 0.029628792824062226
 0.0306227531696316
 ⋮
 0.02665056333747046
 0.02525893438851731
 0.024749347823716496
 0.024158682609370527
 0.023572491294117648
 0.0240760416039372
 0.023230166
 0.023759599761907887
 0.02369791084482974

ros50 = chemicaldf[!, "H2O2 50uM Mean"]
ros50_error = chemicaldf[!, "H2O2 50uM SD"] ./ sqrt.(chemicaldf[!, "H2O2 50uM N"])
ros200 = chemicaldf[!, "H2O2 200uM Mean"]
ros200_error = chemicaldf[!, "H2O2 200uM SD"] ./ sqrt.(chemicaldf[!, "H2O2 200uM N"])

42-element Vector{Float64}:
 0.01579193843196981
 0.015098907892605207
 0.01579878654371588
 0.015408311367567479
 0.01575464714359121
 0.015563615686070896
 0.015839325542753135
 0.01815890968670797
 0.022420341690296358
 0.026401360928455686
 ⋮
 0.01970504059747541
 0.018928610243134615
 0.018308801083333333
 0.017872510075561544
 0.017135324764712814
 0.016622042234321177
 0.01699837223508933
 0.01644948545788867
 0.016236308371972322

plot(ts, ctl, yerr=ctl_error, lab="Control", color=:blue, markerstrokecolor=:blue)
plot!(ts, ros50, yerr=ros50_error, lab="H2O2 50uM", color=:red, markerstrokecolor=:red)
plot!(ts, ros200, yerr=ros200_error, lab="H2O2 200uM", color=:green, markerstrokecolor=:green)
plot!(xlabel="Time (s)", ylabel="CaMKII activity (A.U.)", title="Experiment")

savefig("ros-exp.pdf")

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

Decay rates¶

Fit against an exponential decay model.

decay_model(p, x) = @. p[1] * exp(-x / p[2]) + p[3]

decay_model (generic function with 1 method)

Data from experiments (starts from 120 seconds till the end of the experiment).

xdata = collect(range(0.0, step=5.0, length=18))
ydata_ctl = ctl[25:end]
ydata_50 = ros50[25:end]
ydata_200 = ros200[25:end]
ydata_ctl_sim = sol(stimend:5second:tend, idxs=sys.CaMKAct * 100).u
ydata_01_sim = sol2(stimend:5second:tend, idxs=sys.CaMKAct * 100).u
ydata_05_sim = sol3(stimend:5second:tend, idxs=sys.CaMKAct * 100).u

18-element Vector{Float64}:
 69.01072725422472
 62.165463361824244
 55.3789461717103
 49.953045527169806
 45.50863831977352
 41.796175689293605
 38.64234837936705
 35.92464558959349
 33.55493673973885
 31.4683292740543
 29.61641731445689
 27.96199004614117
 26.476414893863115
 25.13677938179542
 23.924410475946274
 22.824097269700044
 21.823006891318684
 20.91049311643447

Fit data to an exponential decay model

fit_ctl = solve(CurveFitProblem(xdata, ydata_ctl), ExpSumFitAlgorithm(n=1, withconst=true))
fit_50 = solve(CurveFitProblem(xdata, ydata_50), ExpSumFitAlgorithm(n=1, withconst=true))
fit_200 = solve(CurveFitProblem(xdata, ydata_200), ExpSumFitAlgorithm(n=1, withconst=true))
fit_ctl_sim = solve(CurveFitProblem(xdata, ydata_ctl_sim), ExpSumFitAlgorithm(n=1, withconst=true))
fit_01_sim = solve(CurveFitProblem(xdata, ydata_01_sim), ExpSumFitAlgorithm(n=1, withconst=true))
fit_05_sim = solve(CurveFitProblem(xdata, ydata_05_sim), ExpSumFitAlgorithm(n=1, withconst=true))

retcode: Success
alg: ExpSumFitAlgorithm
residuals mean: -5.723816482511918e-15
u: [17.050524250379294, 51.63110628858936, -0.028911808977288095]

Calculate time scales (tau) from fit parameters

tau_exp_ctl = inv(-fit_ctl.u.λ[])
tau_exp_50 = inv(-fit_50.u.λ[])
tau_exp_200 = inv(-fit_200.u.λ[])
tau_sim_ctl = inv(-fit_ctl_sim.u.λ[])
tau_sim_01 = inv(-fit_01_sim.u.λ[])
tau_sim_05 = inv(-fit_05_sim.u.λ[])

println("The time scales for experiments: ")
for (tau, freq) in zip((tau_exp_ctl, tau_exp_50, tau_exp_200), (0, 50, 200))
    println("$freq uM ROS is $(round(tau; digits=2)) seconds.")
end

println("\nThe time scales for simulations: ")
for (tau, freq) in zip((tau_sim_ctl, tau_sim_01, tau_sim_05), (0, 0.1, 0.5))
    println("$freq uM ROS is $(round(tau; digits=2)) seconds.")
end

The time scales for experiments: 
0 uM ROS is 25.82 seconds.
50 uM ROS is 31.79 seconds.
200 uM ROS is 37.46 seconds.

The time scales for simulations: 
0 uM ROS is 20.11 seconds.
0.1 uM ROS is 25.47 seconds.
0.5 uM ROS is 34.59 seconds.

This notebook was generated using Literate.jl.