Figure 6: Calcium oscillation

Figure 6: Calcium oscillation#

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

14

@named sys = make_model()
@unpack Ca_c, GlcConst, kATPCa, kATP = sys
alg = TRBDF2()

TRBDF2(; linsolve = nothing, nlsolve = OrdinaryDiffEqNonlinearSolve.NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}(1//100, 10, 1//5, 1//5, false, true, 0//1), precs = DEFAULT_PRECS, smooth_est = true, extrapolant = linear, controller = PI, step_limiter! = trivial_limiter!,)

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

@named sysosci = make_model(; caceq=cac_wave(amplitude=0.8))
\[\begin{split} \begin{align} \frac{\mathrm{d} \mathrm{NADH}_{m}\left( t \right)}{\mathrm{d}t} &= \frac{J_{DH}\left( t \right) + J_{NADHT}\left( t \right) - J_{O2}\left( t \right)}{V_{MTX}} - kNADHm \mathrm{NADH}_{m}\left( t \right) \\ \frac{\mathrm{d} \mathrm{NADH}_{c}\left( t \right)}{\mathrm{d}t} &= \frac{ - J_{LDH}\left( t \right) - J_{NADHT}\left( t \right) + J_{GPD}\left( t \right)}{V_{I}} - kNADHc \mathrm{NADH}_{c}\left( t \right) \\ \frac{\mathrm{d} \mathrm{Ca}_{m}\left( t \right)}{\mathrm{d}t} &= \frac{F_{M} \left( J_{MCU}\left( t \right) - J_{NCLX}\left( t \right) \right)}{V_{MTX}} \\ \frac{\mathrm{d} \Delta{\Psi}m\left( t \right)}{\mathrm{d}t} &= \frac{ - 2 J_{MCU}\left( t \right) + J_{HR}\left( t \right) - J_{HL}\left( t \right) - J_{HF}\left( t \right) - J_{ANT}\left( t \right)}{C_{MIT}} \\ \frac{\mathrm{d} \mathrm{G3P}\left( t \right)}{\mathrm{d}t} &= \frac{ - J_{GPD}\left( t \right) + 2 J_{GK}\left( t \right)}{V_{I}} - kG3P \mathrm{G3P}\left( t \right) \\ \frac{\mathrm{d} \mathrm{Pyr}\left( t \right)}{\mathrm{d}t} &= \frac{ - J_{LDH}\left( t \right) - J_{PDH}\left( t \right) + J_{GPD}\left( t \right)}{V_{I} + V_{MTX}} - kPyr \mathrm{Pyr}\left( t \right) \\ \frac{\mathrm{d} \mathrm{AEC}\left( t \right)}{\mathrm{d}t} &= \frac{\frac{J_{ANT}\left( t \right) + 2 J_{GPD}\left( t \right) - ATPstiochGK J_{GK}\left( t \right)}{V_{I}} + \left( - kATP - kATPCa \mathrm{Ca}_{c}\left( t \right) \right) \mathrm{ATP}_{c}\left( t \right)}{{\Sigma}Ac} \\ \frac{\mathrm{d} \mathrm{x2}\left( t \right)}{\mathrm{d}t} &= - \mathrm{tipside}\left( t \right) + \mathrm{tiptip}\left( t \right) \\ \frac{\mathrm{d} \mathrm{x3}\left( t \right)}{\mathrm{d}t} &= \mathrm{tipside}\left( t \right) \end{align} \end{split}\]
tend = 4000.0
ts = range(tend-480, tend; length=201)
prob = ODEProblem(sysosci, [], tend, [GlcConst => 10mM])
sol = solve(prob, alg, saveat=ts)

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.10663925370712568, 0.0031096626353270935, 0.0008340480714225096, 0.12216920919929514, 0.005308719248177244, 0.034742727964047364, 0.9825657172782265, 0.21294081412927773, 0.03089426697951148]
 [0.10715380690566202, 0.003120781524156824, 0.0008280597985221397, 0.12422917394683575, 0.00531389221966499, 0.034922944276474986, 0.9835691005305924, 0.21322745963350612, 0.03097965784346304]
 [0.10777112641020663, 0.003133981571092565, 0.0008189726388351657, 0.12658002124297452, 0.0053196450682061235, 0.03514482204388877, 0.9845936021690224, 0.21343764547694055, 0.031044026884199515]
 [0.1084841088329694, 0.003149075283636527, 0.0008065875422923966, 0.1290805252802413, 0.005325554502885413, 0.03540894195115609, 0.9855315995606376, 0.2135710990034269, 0.03108675537909771]
 [0.10926966264649839, 0.003165587631750429, 0.0007909550198356732, 0.13161741967180737, 0.005331357366361329, 0.035710880400487754, 0.9863514370952637, 0.21363167896029212, 0.031108764416959138]
 [0.11010337105684423, 0.003183027284557047, 0.0007721829202800557, 0.13409908424749947, 0.005336865620302413, 0.036045460167016444, 0.9870453365627672, 0.21362363851519334, 0.031111190510053576]
 [0.11093646933120761, 0.0032003987945438965, 0.0007507205065172073, 0.13636930602132252, 0.005341770829424758, 0.03639998718967605, 0.9876059942472781, 0.21355591086544384, 0.031096488078505832]
 [0.11169364023857316, 0.0032161259043552255, 0.0007277245273543962, 0.13824320977869956, 0.005345730666318264, 0.03674887600576968, 0.9880388672661268, 0.21344526348567752, 0.031069359198613124]
 [0.11237147364622882, 0.003230090889240446, 0.0007036703125888848, 0.13976743070243916, 0.0053488533746290385, 0.037084162659176315, 0.9883599524951295, 0.2132973053388606, 0.031031404162040883]
 [0.1129673379833741, 0.0032422367999558254, 0.000679275906446586, 0.1410081740493529, 0.005351303941245797, 0.03739643571097978, 0.9885954289224343, 0.2131195020160714, 0.030984759420997637]
 ⋮
 [0.10600143580491707, 0.003092839833581943, 0.0007583225542293416, 0.11684793416555404, 0.005294343523070257, 0.035027469305702184, 0.979102000922447, 0.20862976081454587, 0.029696591933692583]
 [0.10568204602047501, 0.0030866450993007454, 0.000783934167304157, 0.11632462356242258, 0.00529290272217886, 0.03473407714307588, 0.978726695840973, 0.209270651866145, 0.02986716872245448]
 [0.1055106671379244, 0.003083434726534486, 0.0008033508041541498, 0.1161785242850786, 0.005292553902015528, 0.034554920352817944, 0.9786695545159343, 0.209911105220551, 0.03003999595098704]
 [0.10546371235861747, 0.0030828165251904054, 0.0008177734174494727, 0.11641216037967349, 0.0052932332383926715, 0.03445551219227657, 0.9788780368232399, 0.21053008212992633, 0.030209224713816575]
 [0.1055183950146881, 0.0030844063209235746, 0.0008279445838303766, 0.11696597111103389, 0.005294806091239262, 0.0344183438356316, 0.9793115162099723, 0.21111512913771727, 0.030371263270650307]
 [0.10566740905906201, 0.0030880077661655132, 0.0008341925017411616, 0.11782577959474122, 0.0052972161582826695, 0.0344368555167276, 0.9799469432063266, 0.211660401491637, 0.030524358711425708]
 [0.10590347889549122, 0.003093417050902516, 0.0008371048235023548, 0.1189788133168179, 0.0053003703126559935, 0.03450034877178243, 0.9807277602826139, 0.21215432570022055, 0.0306649313358123]
 [0.10622443754490846, 0.0031005781373967148, 0.000837095943658422, 0.1204226898800985, 0.005304208489228372, 0.03460164161517213, 0.9816089019054227, 0.21258592795711684, 0.030789467666632982]
 [0.10663684114372977, 0.0031096166990679888, 0.0008341759394805238, 0.12217574638737118, 0.005308749015833915, 0.034740951582228015, 0.9825778969584819, 0.21295014886260577, 0.030896247834105738]

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)
_images/0e0d3270c1aac9a241654cf4fdaebd0bdba059cba840193474337249f970e5fe.png

Export figure

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

Python: None

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

prob2 = ODEProblem(sysosci, [], tend, [GlcConst => 10mM, kATPCa=>10Hz/mM, kATP=>0.055Hz])
sol2 = solve(prob2, alg, saveat=ts)
plot_fig5(sol2)
_images/2923a2abe4e0c4d6b0fbb3569c59d043fc22ff18efb85987f135ba56efb41802.png

kATPCa : 90 -> 0.1

prob4 = ODEProblem(sysosci, [], tend, [GlcConst => 10mM, kATPCa=>0.1Hz/mM, kATP=>0.06Hz])
sol4 = solve(prob4, alg, saveat=ts)
plot_fig5(sol4)
_images/bc63931936f27f10335fff25b35ec5676eb17d12f5906825664128eeeb486c88.png

This notebook was generated using Literate.jl.