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

14

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

Build system: 150.033447 seconds (107.16 M allocations: 5.196 GiB, 2.95% gc time, 99.96% 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.751828 seconds (2.36 M allocations: 118.339 MiB, 98.29% 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: 33.388448 seconds (61.01 M allocations: 2.985 GiB, 4.94% gc time, 99.84% compilation time: 10% of which was recompilation)
Solve oscillating problem: 15.522339 seconds (18.73 M allocations: 924.693 MiB, 2.02% gc time, 99.80% 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.030897730017074507, 0.21295121162387848, 0.9825648821436483, 0.03474278577058455, 0.00530871945675025, 0.12216997896651598, 0.0008339595522119663, 0.003109667789480741, 0.10663953948545156]
 [0.030983154491593518, 0.21323782989569706, 0.9835708365128284, 0.0349228832994455, 0.0053138973228006, 0.12423060526308118, 0.0008279955103211213, 0.003120779373961228, 0.10715364768942368]
 [0.031047568835679572, 0.21344808428341058, 0.9845945677261049, 0.03514464683575063, 0.005319641182627445, 0.12657733610403582, 0.0008189240346874688, 0.0031339653102348468, 0.10777025781686739]
 [0.031090334166617398, 0.2135816702238589, 0.9855286453767145, 0.03540857385350167, 0.005325531757812094, 0.12907086057489797, 0.0008065454086560407, 0.0031490368242044903, 0.10848228953810621]
 [0.03111239096604965, 0.21364239850794256, 0.9863483679852519, 0.03571025692671435, 0.005331332896040528, 0.1316070881479704, 0.0007909217989272775, 0.003165536169907016, 0.10926731147002895]
 [0.031114833348483372, 0.2136343806845818, 0.9870440722142324, 0.0360447861570321, 0.005336853860804066, 0.13409460046648902, 0.0007721531241019146, 0.0031829921887963544, 0.11010187849519218]
 [0.0311001245409317, 0.21356659358257854, 0.9876062748686222, 0.036399424982840774, 0.005341772414467536, 0.13637107662953182, 0.0007506934818901028, 0.0032003923256828585, 0.11093630947539844]
 [0.031072961987527084, 0.2134557985657959, 0.9880392462969871, 0.03674854571293962, 0.005345733102832202, 0.138244765282835, 0.0007276980850537018, 0.0032161281015491484, 0.11169369708993133]
 [0.031034945019485224, 0.21330759976931937, 0.9883596825603593, 0.03708411784102084, 0.005348850978364386, 0.13976661723616385, 0.0007036380984835844, 0.0032300903304615745, 0.11237136494094102]
 [0.03098828949961417, 0.21312973889771492, 0.9885952818978386, 0.03739639314762703, 0.005351300502996882, 0.14100708904024895, 0.0006792569659526862, 0.0032422298082991543, 0.11296701529843278]
 ⋮
 [0.02969694577772109, 0.20863287482886728, 0.9791092412806476, 0.035032595968497865, 0.005294406293249448, 0.11687290355364599, 0.0007581414384417866, 0.003093089128992541, 0.10601039901911924]
 [0.029867097014710065, 0.20927214578659717, 0.9787246938722692, 0.034734932820572245, 0.005292885635492483, 0.11632000029765233, 0.0007839010217544078, 0.003086591978501745, 0.10568189648263764]
 [0.030038975656994284, 0.20990910343250405, 0.9786628438025768, 0.03454974606740217, 0.005292533590407104, 0.11617250538335382, 0.0008035082044554853, 0.003083341138772003, 0.10550587811057154]
 [0.030209185114656263, 0.21053143533206847, 0.9788882106015029, 0.034456111966831036, 0.00529327048736656, 0.11642372838515619, 0.0008176907831598212, 0.0030829012238427268, 0.10546600770658707]
 [0.03037193427487195, 0.21111885793639448, 0.9793244798961159, 0.03442231770695518, 0.005294855136153329, 0.11698117930893326, 0.0008277404583113251, 0.0030845419645666745, 0.10552386542618826]
 [0.03052510144546169, 0.21166431101264296, 0.9799496910742072, 0.034440274745273135, 0.005297236211846659, 0.11783243400722379, 0.0008340330693120032, 0.003088089046995444, 0.10567140634561281]
 [0.030665554906546624, 0.2121578233137819, 0.9807260454046959, 0.03450234977574106, 0.005300367282516573, 0.11897825246214001, 0.0008369979257212281, 0.003093443491281408, 0.10590518034293435]
 [0.030789551034468808, 0.21258781505169838, 0.9816049294128835, 0.0346027973052286, 0.00530420219290503, 0.12042295534768986, 0.0008370132926381447, 0.0031006035598615815, 0.1062260504172179]
 [0.03089581632986506, 0.2129505330168972, 0.9825753611180086, 0.03474175674585007, 0.0053087374073961715, 0.12217138247762521, 0.0008340906929471155, 0.0031096216566007355, 0.1066368207234749]

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.007839 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.006207 seconds (1.01 k allocations: 101.477 KiB)

This notebook was generated using Literate.jl.