Figure 6: Calcium oscillation - Mitochondrial dynamics

using OrdinaryDiffEq
using OrdinaryDiffEqSDIRK
using ModelingToolkit
using MitochondrialDynamics
using MitochondrialDynamics: second, μM, mV, mM, Hz, minute, hil
import PythonPlot as plt
plt.matplotlib.rcParams["font.size"] = 14

┌ Info: CondaPkg: Waiting for lock to be freed. You may delete this file if no other process is resolving.
└   lock_file = "/home/github/actions-runner-1/_work/mitodyn-ode/mitodyn-ode/.CondaPkg/lock"

14

@time "Build system" @named sys = make_model()
@unpack Ca_c, Glc, kATPCa, kATP = sys
alg = TRBDF2()

Build system: 87.091377 seconds (108.31 M allocations: 5.253 GiB, 3.73% gc time, 99.51% compilation time: 4% of which was recompilation)

TRBDF2(; linsolve = nothing, nlsolve = OrdinaryDiffEqNonlinearSolve.NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Nothing}(1//100, 10, 1//5, 1//5, false, true, nothing), smooth_est = true, predictor = Linear, step_limiter! = trivial_limiter!, autodiff = ADTypes.AutoForwardDiff(), concrete_jac = nothing,)

Calcium oscillation function

function cac_wave(; ca_base = 0.09μM, ca_act = 0.25μM, n=4, katp=25, amplitude=0.5, period=2minute)
    @variables t Ca_c(t) ATP_c(t) ADP_c(t)
    @parameters (RestingCa=ca_base, ActivatedCa=ca_act, NCac=n, KatpCac=katp)
    x = 5 * ((t / period) % 1.0) ## An oscillating function
    w = (x * exp(1 - x))^4  ## Scale from 0 to 1
    caceq = Ca_c ~ RestingCa + ActivatedCa * hil(ATP_c, KatpCac * ADP_c, NCac) * (1 + amplitude * (2w-1))
    return caceq
end

@time "Build oscillating system" @named sysosci = make_model(; caceq=cac_wave(amplitude=0.8))

Build oscillating system: 1.044500 seconds (2.36 M allocations: 118.360 MiB, 4.31% gc time, 98.01% compilation time: 2% of which was recompilation)

Model sysosci:
Equations (9):
  9 standard: see equations(sysosci)
Unknowns (9): see unknowns(sysosci)
  x3(t)
  x2(t)
  AEC(t)
  Pyr(t)
  ⋮
Parameters (65): see parameters(sysosci)
  RestingCa
  ActivatedCa
  NCac
  KatpCac
  ⋮
Observed (24): see observed(sysosci)

tend = 4000.0
ts = range(tend-480, tend; length=201)
@time "Build oscillating problem" prob = ODEProblem(sysosci, [], tend, [Glc => 10mM])
@time "Solve oscillating problem" sol = solve(prob, alg, saveat=ts)

┌ Warning: `SciMLBase.ODEProblem(sys, u0, tspan, p; kw...)` is deprecated. Use
│ `SciMLBase.ODEProblem(sys, merge(if isempty(u0)
│     Dict()
│ else
│     Dict(unknowns(sys) .=> u0)
│ end, Dict(p)), tspan)` instead.
└ @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/PHXso/src/deprecations.jl:53
Build oscillating problem: 32.091410 seconds (61.01 M allocations: 2.986 GiB, 3.37% gc time, 99.85% compilation time: 11% of which was recompilation)
Solve oscillating problem: 15.205594 seconds (18.77 M allocations: 925.377 MiB, 5.40% gc time, 99.81% compilation time)

retcode: Success
Interpolation: 1st order linear
t: 201-element Vector{Float64}:
 3520.0
 3522.4
 3524.8
 3527.2
 3529.6
 3532.0
 3534.4
 3536.8
 3539.2
 3541.6
    ⋮
 3980.8
 3983.2
 3985.6
 3988.0
 3990.4
 3992.8
 3995.2
 3997.6
 4000.0
u: 201-element Vector{Vector{Float64}}:
 [0.03089773001707452, 0.21295121162387842, 0.9825648821436481, 0.03474278577058456, 0.00530871945675025, 0.12216997896651598, 0.0008339595522119662, 0.0031096677894807425, 0.10663953948545159]
 [0.030983154491593532, 0.213237829895697, 0.9835708365128282, 0.0349228832994455, 0.005313897322800599, 0.12423060526308116, 0.0008279955103211212, 0.0031207793739612285, 0.10715364768942369]
 [0.031047568835679586, 0.21344808428341053, 0.9845945677261048, 0.03514464683575063, 0.0053196411826274435, 0.1265773361040358, 0.0008189240346874687, 0.003133965310234847, 0.10777025781686736]
 [0.03109033416661741, 0.21358167022385885, 0.9855286453767144, 0.03540857385350168, 0.005325531757812093, 0.12907086057489794, 0.0008065454086560406, 0.0031490368242044903, 0.1084822895381062]
 [0.03111239096604966, 0.21364239850794253, 0.9863483679852518, 0.03571025692671435, 0.005331332896040528, 0.13160708814797045, 0.0007909217989272773, 0.003165536169907016, 0.10926731147002894]
 [0.031114833348483383, 0.21363438068458174, 0.9870440722142323, 0.0360447861570321, 0.005336853860804066, 0.13409460046648908, 0.0007721531241019144, 0.003182992188796353, 0.11010187849519218]
 [0.03110012454093171, 0.21356659358257848, 0.9876062748686222, 0.03639942498284079, 0.005341772414467535, 0.1363710766295318, 0.0007506934818901026, 0.0032003923256828564, 0.11093630947539844]
 [0.031072961987527094, 0.21345579856579586, 0.988039246296987, 0.03674854571293963, 0.005345733102832202, 0.1382447652828349, 0.0007276980850537012, 0.003216128101549147, 0.11169369708993135]
 [0.031034945019485234, 0.21330759976931932, 0.9883596825603594, 0.03708411784102086, 0.005348850978364386, 0.1397666172361638, 0.0007036380984835838, 0.0032300903304615745, 0.11237136494094103]
 [0.030988289499614182, 0.21312973889771483, 0.9885952818978387, 0.03739639314762704, 0.005351300502996882, 0.14100708904024892, 0.0006792569659526855, 0.0032422298082991547, 0.11296701529843275]
 ⋮
 [0.02969694577772109, 0.20863287482886722, 0.9791092412806476, 0.03503259596849786, 0.005294406293249448, 0.11687290355364606, 0.0007581414384417879, 0.0030930891289925426, 0.10601039901911928]
 [0.029867097014710065, 0.2092721457865971, 0.9787246938722693, 0.03473493282057225, 0.005292885635492485, 0.1163200002976525, 0.0007839010217544088, 0.003086591978501747, 0.1056818964826377]
 [0.030038975656994284, 0.209909103432504, 0.9786628438025771, 0.034549746067402194, 0.005292533590407104, 0.1161725053833539, 0.0008035082044554862, 0.0030833411387720046, 0.1055058781105716]
 [0.030209185114656263, 0.21053143533206845, 0.9788882106015029, 0.034456111966831064, 0.005293270487366561, 0.1164237283851563, 0.0008176907831598223, 0.0030829012238427276, 0.10546600770658716]
 [0.03037193427487195, 0.21111885793639443, 0.9793244798961158, 0.03442231770695521, 0.00529485513615333, 0.11698117930893334, 0.0008277404583113258, 0.0030845419645666745, 0.1055238654261883]
 [0.03052510144546169, 0.2116643110126429, 0.979949691074207, 0.03444027474527315, 0.0052972362118466594, 0.11783243400722382, 0.0008340330693120037, 0.0030880890469954436, 0.10567140634561283]
 [0.030665554906546624, 0.21215782331378186, 0.9807260454046958, 0.03450234977574106, 0.005300367282516571, 0.11897825246214001, 0.0008369979257212286, 0.003093443491281408, 0.1059051803429344]
 [0.030789551034468808, 0.21258781505169833, 0.9816049294128835, 0.03460279730522858, 0.005304202192905029, 0.12042295534768982, 0.0008370132926381448, 0.0031006035598615794, 0.1062260504172179]
 [0.03089581632986506, 0.21295053301689715, 0.9825753611180086, 0.034741756745850054, 0.0053087374073961715, 0.12217138247762521, 0.0008340906929471154, 0.003109621656600734, 0.10663682072347488]

function plot_fig5(sol, figsize=(10, 10))
    ts = sol.t
    tsm = ts ./ 60
    @unpack Ca_c, Ca_m, ATP_c, ADP_c, ΔΨm, degavg, J_ANT, J_HL = sys
    fig, ax = plt.subplots(5, 1; figsize)

    ax[0].plot(tsm, sol[Ca_c * 1000], label="Cyto")
    ax[0].plot(tsm, sol[Ca_m * 1000], label="Mito")
    ax[0].set_title("a", loc="left")
    ax[0].set_ylabel("Calcium (μM)", fontsize=12)
    ax[0].legend(loc="center left")

    ax[1].plot(tsm, sol[ATP_c / ADP_c])
    ax[1].set_title("b", loc="left")
    ax[1].set_ylabel("ATP:ADP (ratio)", fontsize=12)

    ax[2].plot(tsm, sol[ΔΨm * 1000])
    ax[2].set_title("c", loc="left")
    ax[2].set_ylabel("ΔΨm (mV)", fontsize=12)

    ax[3].plot(tsm, sol[degavg], label="Average node degree")
    ax[3].set_title("d", loc="left")
    ax[3].set_ylabel("a.u.")
    ax[3].legend(loc="center left")

    ax[4].plot(tsm, sol[J_ANT], label="ATP export")
    ax[4].plot(tsm, sol[J_HL], label="H leak")
    ax[4].set_title("e", loc="left")
    ax[4].set_ylabel("Rate (mM/s)")
    ax[4].set(xlabel="Time (minute)")
    ax[4].legend(loc="center left")

    for i in 0:4
        ax[i].grid()
        ax[i].set_xlim(tsm[begin], tsm[end])
    end

    plt.tight_layout()
    return fig
end

plot_fig5 (generic function with 2 methods)

fig5 = plot_fig5(sol)

Python: <Figure size 1000x1000 with 5 Axes>

Export figure

exportTIF(fig5, "Fig6-ca-oscillation.tif")

Python: None

Tuning ca-dependent ATP consumption rate (kATPCa) kATPCa : 90 -> 10

prob2 = ODEProblem(sysosci, [], tend, [Glc => 10mM, kATPCa=>10Hz/mM, kATP=>0.055Hz])
@time "Solve tuned problem 2" sol2 = solve(prob2, alg, saveat=ts)
plot_fig5(sol2)

┌ Warning: `SciMLBase.ODEProblem(sys, u0, tspan, p; kw...)` is deprecated. Use
│ `SciMLBase.ODEProblem(sys, merge(if isempty(u0)
│     Dict()
│ else
│     Dict(unknowns(sys) .=> u0)
│ end, Dict(p)), tspan)` instead.
└ @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/PHXso/src/deprecations.jl:53
Solve tuned problem 2: 0.004113 seconds (1.01 k allocations: 101.523 KiB)

kATPCa : 90 -> 0.1

prob4 = ODEProblem(sysosci, [], tend, [Glc => 10mM, kATPCa=>0.1Hz/mM, kATP=>0.06Hz])
@time "Solve tuned problem 4" sol4 = solve(prob4, alg, saveat=ts)
plot_fig5(sol4)

┌ Warning: `SciMLBase.ODEProblem(sys, u0, tspan, p; kw...)` is deprecated. Use
│ `SciMLBase.ODEProblem(sys, merge(if isempty(u0)
│     Dict()
│ else
│     Dict(unknowns(sys) .=> u0)
│ end, Dict(p)), tspan)` instead.
└ @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/PHXso/src/deprecations.jl:53
Solve tuned problem 4: 0.003758 seconds (1.01 k allocations: 101.477 KiB)

This notebook was generated using Literate.jl.