Figure 6: Calcium oscillation

Figure 6: Calcium oscillation#

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

┌ Warning: Module SymbolicsPreallocationToolsExt with build ID ffffffff-ffff-ffff-8a55-7cdb31025361 is missing from the cache.
│ This may mean SymbolicsPreallocationToolsExt [d479e226-fb54-5ebe-a75e-a7af7f39127f] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:2018
┌ Error: Error during loading of extension SymbolicsPreallocationToolsExt of Symbolics, use `Base.retry_load_extensions()` to retry.
│   exception =
│    1-element ExceptionStack:
│    Declaring __precompile__(false) is not allowed in files that are being precompiled.
│    Stacktrace:
│      [1] _require(pkg::Base.PkgId, env::Nothing)
│        @ Base ./loading.jl:2022
│      [2] __require_prelocked(uuidkey::Base.PkgId, env::Nothing)
│        @ Base ./loading.jl:1882
│      [3] #invoke_in_world#3
│        @ ./essentials.jl:926 [inlined]
│      [4] invoke_in_world
│        @ ./essentials.jl:923 [inlined]
│      [5] _require_prelocked
│        @ ./loading.jl:1873 [inlined]
│      [6] _require_prelocked
│        @ ./loading.jl:1872 [inlined]
│      [7] run_extension_callbacks(extid::Base.ExtensionId)
│        @ Base ./loading.jl:1365
│      [8] run_extension_callbacks(pkgid::Base.PkgId)
│        @ Base ./loading.jl:1400
│      [9] run_package_callbacks(modkey::Base.PkgId)
│        @ Base ./loading.jl:1224
│     [10] __require_prelocked(uuidkey::Base.PkgId, env::String)
│        @ Base ./loading.jl:1889
│     [11] #invoke_in_world#3
│        @ ./essentials.jl:926 [inlined]
│     [12] invoke_in_world
│        @ ./essentials.jl:923 [inlined]
│     [13] _require_prelocked(uuidkey::Base.PkgId, env::String)
│        @ Base ./loading.jl:1873
│     [14] macro expansion
│        @ ./loading.jl:1860 [inlined]
│     [15] macro expansion
│        @ ./lock.jl:267 [inlined]
│     [16] __require(into::Module, mod::Symbol)
│        @ Base ./loading.jl:1823
│     [17] #invoke_in_world#3
│        @ ./essentials.jl:926 [inlined]
│     [18] invoke_in_world
│        @ ./essentials.jl:923 [inlined]
│     [19] require(into::Module, mod::Symbol)
│        @ Base ./loading.jl:1816
│     [20] include
│        @ ./Base.jl:495 [inlined]
│     [21] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt128}}, source::String)
│        @ Base ./loading.jl:2292
│     [22] top-level scope
│        @ stdin:4
│     [23] eval
│        @ ./boot.jl:385 [inlined]
│     [24] include_string(mapexpr::typeof(identity), mod::Module, code::String, filename::String)
│        @ Base ./loading.jl:2146
│     [25] include_string
│        @ ./loading.jl:2156 [inlined]
│     [26] exec_options(opts::Base.JLOptions)
│        @ Base ./client.jl:321
│     [27] _start()
│        @ Base ./client.jl:557
└ @ Base loading.jl:1371

14

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

TRBDF2(; 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, step_limiter! = trivial_limiter!, autodiff = ADTypes.AutoForwardDiff(),)

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

equations(sysosci)

observed(sysosci)
\[\begin{split} \begin{align} \mathtt{x1}\left( t \right) &= 1 - 3 \mathtt{x3}\left( t \right) - 2 \mathtt{x2}\left( t \right) \\ \mathtt{ATP\_c}\left( t \right) &= \left( \frac{1 - \sqrt{1 + 4 \left( -1 + 4 \mathtt{kEqAK} \right) \mathtt{AEC}\left( t \right) \left( 1 - \mathtt{AEC}\left( t \right) \right)}}{-2 + 8 \mathtt{kEqAK}} + \mathtt{AEC}\left( t \right) \right) \mathtt{{\Sigma}Ac} \\ \mathtt{ADP\_c}\left( t \right) &= \frac{2 \left( -1 + \sqrt{1 + 4 \left( -1 + 4 \mathtt{kEqAK} \right) \mathtt{AEC}\left( t \right) \left( 1 - \mathtt{AEC}\left( t \right) \right)} \right) \mathtt{{\Sigma}Ac}}{-2 + 8 \mathtt{kEqAK}} \\ \mathtt{J\_HL}\left( t \right) &= \mathtt{pHleak} \mathtt{rHL} e^{\mathtt{kvHleak} \mathtt{\Delta{\Psi}m}\left( t \right)} \\ \mathtt{NAD\_c}\left( t \right) &= - \mathtt{NADH\_c}\left( t \right) + \mathtt{{\Sigma}n\_c} \\ \mathtt{NAD\_m}\left( t \right) &= - \mathtt{NADH\_m}\left( t \right) + \mathtt{{\Sigma}n\_m} \\ \mathtt{J\_HR}\left( t \right) &= \frac{\left( 1 + \mathtt{KaETC} \mathtt{\Delta{\Psi}m}\left( t \right) \right) \mathtt{VmaxETC} \mathtt{rETC} \mathtt{NADH\_m}\left( t \right)}{\left( 1 + \mathtt{KbETC} \mathtt{\Delta{\Psi}m}\left( t \right) \right) \left( \mathtt{KnadhETC} + \mathtt{NADH\_m}\left( t \right) \right)} \\ \mathtt{degavg}\left( t \right) &= \frac{3 \mathtt{x3}\left( t \right) + \mathtt{x1}\left( t \right) + 2 \mathtt{x2}\left( t \right)}{\mathtt{x3}\left( t \right) + \mathtt{x1}\left( t \right) + \mathtt{x2}\left( t \right)} \\ \mathtt{J\_GK}\left( t \right) &= \frac{\mathtt{VmaxGK} \mathtt{ATP\_c}\left( t \right) pow\left( \mathtt{Glc}\left( t \right), \mathtt{nGK} \right)}{\left( \mathtt{KatpGK} + \mathtt{ATP\_c}\left( t \right) \right) \left( pow\left( \mathtt{KglcGK}, \mathtt{nGK} \right) + pow\left( \mathtt{Glc}\left( t \right), \mathtt{nGK} \right) \right)} \\ \mathtt{Ca\_c}\left( t \right) &= \mathtt{RestingCa} + \frac{\mathtt{ActivatedCa} \left( 1 + 0.8 \left( -1 + 1250 \left( rem\left( 0.0083333 t, 1 \right) \right)^{4} \left( e^{1 - 5 rem\left( 0.0083333 t, 1 \right)} \right)^{4} \right) \right) pow\left( \mathtt{ATP\_c}\left( t \right), \mathtt{NCac} \right)}{pow\left( \mathtt{KatpCac} \mathtt{ADP\_c}\left( t \right), \mathtt{NCac} \right) + pow\left( \mathtt{ATP\_c}\left( t \right), \mathtt{NCac} \right)} \\ \mathtt{J\_HF}\left( t \right) &= \frac{\mathtt{VmaxF1} \mathtt{rF1} pow\left( \mathtt{FmgadpF1} \mathtt{ADP\_c}\left( t \right), \mathtt{nadpF1} \right) \left( 1 - e^{\frac{ - \mathtt{Ca\_m}\left( t \right)}{\mathtt{KcaF1}}} \right) pow\left( \mathtt{\Delta{\Psi}m}\left( t \right), \mathtt{nvF1} \right)}{\left( pow\left( \mathtt{KvF1}, \mathtt{nvF1} \right) + pow\left( \mathtt{\Delta{\Psi}m}\left( t \right), \mathtt{nvF1} \right) \right) \left( pow\left( \mathtt{KadpF1}, \mathtt{nadpF1} \right) + pow\left( \mathtt{FmgadpF1} \mathtt{ADP\_c}\left( t \right), \mathtt{nadpF1} \right) \right)} \\ \mathtt{AMP\_c}\left( t \right) &= - \mathtt{ATP\_c}\left( t \right) - \mathtt{ADP\_c}\left( t \right) + \mathtt{{\Sigma}Ac} \\ \mathtt{J\_LDH}\left( t \right) &= \frac{\mathtt{VmaxLDH} \mathtt{Pyr}\left( t \right) \mathtt{NADH\_c}\left( t \right)}{\left( \mathtt{NADH\_c}\left( t \right) + \mathtt{KnadhLDH} \mathtt{NAD\_c}\left( t \right) \right) \left( \mathtt{KpyrLDH} + \mathtt{Pyr}\left( t \right) \right)} \\ \mathtt{J\_GPD}\left( t \right) &= \frac{\mathtt{VmaxGPD} \mathtt{NAD\_c}\left( t \right) \mathtt{G3P}\left( t \right) \mathtt{ADP\_c}\left( t \right)}{\left( \mathtt{KadpGPD} + \mathtt{ADP\_c}\left( t \right) \right) \left( \mathtt{Kg3pGPD} + \mathtt{G3P}\left( t \right) \right) \left( \mathtt{NAD\_c}\left( t \right) + \mathtt{KnadGPD} \mathtt{NADH\_c}\left( t \right) \right)} \\ \mathtt{J\_FFA}\left( t \right) &= \mathtt{kFFA} \mathtt{NAD\_m}\left( t \right) \\ \mathtt{J\_PDH}\left( t \right) &= \frac{\mathtt{VmaxPDH} \mathtt{rPDH} \mathtt{NAD\_m}\left( t \right) \mathtt{Pyr}\left( t \right)}{\left( \mathtt{NAD\_m}\left( t \right) + \mathtt{KnadPDH} \mathtt{NADH\_m}\left( t \right) \right) \left( \mathtt{KpyrPDH} + \mathtt{Pyr}\left( t \right) \right) \left( 1 + \left( 1 + \left( \frac{\mathtt{KcaPDH}}{\mathtt{KcaPDH} + \mathtt{Ca\_m}\left( t \right)} \right)^{2} \mathtt{U1PDH} \right) \mathtt{U2PDH} \right)} \\ \mathtt{J\_NADHT}\left( t \right) &= \frac{\mathtt{VmaxNADHT} \mathtt{NAD\_m}\left( t \right) \mathtt{NADH\_c}\left( t \right)}{\left( \mathtt{NADH\_c}\left( t \right) + \mathtt{Ktn\_c} \mathtt{NAD\_c}\left( t \right) \right) \left( \mathtt{NAD\_m}\left( t \right) + \mathtt{Ktn\_m} \mathtt{NADH\_m}\left( t \right) \right)} \\ \mathtt{J\_O2}\left( t \right) &= \frac{1}{10} \mathtt{J\_HR}\left( t \right) \\ \mathtt{J\_MCU}\left( t \right) &= \frac{74.872 \mathtt{PcaMCU} \left( - 0.2 \mathtt{Ca\_m}\left( t \right) + 0.341 \left( 1 + expm1\left( 74.872 \mathtt{\Delta{\Psi}m}\left( t \right) \right) \right) \mathtt{Ca\_c}\left( t \right) \right) \mathtt{\Delta{\Psi}m}\left( t \right)}{expm1\left( 74.872 \mathtt{\Delta{\Psi}m}\left( t \right) \right)} \\ \mathtt{J\_NCLX}\left( t \right) &= \frac{\mathtt{VmaxNCLX} \left( \frac{\left( \frac{\mathtt{Na\_c}}{\mathtt{KnaNCLX}} \right)^{2} \mathtt{Ca\_m}\left( t \right)}{\mathtt{KcaNCLX}} + \frac{ - \left( \frac{\mathtt{Na\_m}}{\mathtt{KnaNCLX}} \right)^{2} \mathtt{Ca\_c}\left( t \right)}{\mathtt{KcaNCLX}} \right)}{1 + \frac{\mathtt{Ca\_m}\left( t \right)}{\mathtt{KcaNCLX}} + \frac{\left( \frac{\mathtt{Na\_c}}{\mathtt{KnaNCLX}} \right)^{2} \mathtt{Ca\_m}\left( t \right)}{\mathtt{KcaNCLX}} + \frac{\mathtt{Ca\_c}\left( t \right)}{\mathtt{KcaNCLX}} + \frac{\left( \frac{\mathtt{Na\_m}}{\mathtt{KnaNCLX}} \right)^{2} \mathtt{Ca\_c}\left( t \right)}{\mathtt{KcaNCLX}} + \left( \frac{\mathtt{Na\_c}}{\mathtt{KnaNCLX}} \right)^{2} + \left( \frac{\mathtt{Na\_m}}{\mathtt{KnaNCLX}} \right)^{2}} \\ \mathtt{J\_ANT}\left( t \right) &= \frac{1}{3} \mathtt{J\_HF}\left( t \right) \\ \mathtt{J\_DH}\left( t \right) &= \mathtt{J\_FFA}\left( t \right) + 4.6 \mathtt{J\_PDH}\left( t \right) \\ \mathtt{tipside}\left( t \right) &= \frac{\mathtt{kfuse2} \mathtt{J\_ANT}\left( t \right) \mathtt{x1}\left( t \right) \mathtt{x2}\left( t \right)}{\mathtt{J\_HL}\left( t \right)} - \mathtt{kfiss2} \mathtt{x3}\left( t \right) \\ \mathtt{tiptip}\left( t \right) &= \frac{\left( \mathtt{x1}\left( t \right) \right)^{2} \mathtt{kfuse1} \mathtt{J\_ANT}\left( t \right)}{\mathtt{J\_HL}\left( t \right)} - \mathtt{kfiss1} \mathtt{x2}\left( t \right) \end{align} \end{split}\]
tend = 4000.0
ts = range(tend-480, tend; length=201)
prob = ODEProblem(sysosci, [], tend, [Glc => 10mM])
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.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/hlGut/src/deprecations.jl:45

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.07298768624306846, 0.24897140514600974, 0.8383275916263327, 0.031581882998430205, 0.005079515343908573, 0.09807633528770028, 0.00021267386427249312, 0.0028124343488794994, 0.09905459455316731]
 [0.07298765519282241, 0.2489713650818988, 0.8383273436364553, 0.031581860740286294, 0.005079514977486805, 0.09807630785119424, 0.00021267452146499945, 0.002812433784045541, 0.09905457662206517]
 [0.07298762486822354, 0.24897132715525125, 0.8383272196614615, 0.031581843441071615, 0.005079514778174247, 0.09807628978868912, 0.00021267490775616927, 0.0028124333970135024, 0.09905456148956045]
 [0.07298759525189746, 0.24897129138259846, 0.8383272137796625, 0.03158183093548212, 0.005079514738094954, 0.09807628067993324, 0.00021267503526340743, 0.0028124331804524, 0.09905454908564722]
 [0.07298756632646984, 0.24897125778047188, 0.8383273200693698, 0.03158182305821381, 0.005079514849372986, 0.09807628010467498, 0.00021267491610411863, 0.0028124331270312497, 0.09905453934031959]
 [0.0729875380745663, 0.24897122636540278, 0.8383275326088951, 0.03158181964396263, 0.005079515104132401, 0.09807628764266269, 0.00021267456239570769, 0.0028124332294190684, 0.09905453218357171]
 [0.07298751047881247, 0.24897119715392263, 0.8383278454765498, 0.03158182052742458, 0.005079515494497256, 0.09807630287364472, 0.00021267398625557932, 0.002812433480284872, 0.0990545275453977]
 [0.07298748352183401, 0.24897117016256273, 0.8383282527506458, 0.03158182554329563, 0.00507951601259161, 0.09807632537736945, 0.00021267319980113842, 0.0028124338722976764, 0.09905452535579168]
 [0.07298745718625653, 0.24897114540785453, 0.8383287485094943, 0.03158183452627176, 0.00507951665053952, 0.09807635473358522, 0.00021267221514978963, 0.0028124343981264983, 0.09905452554474774]
 [0.07298743145470567, 0.24897112290632942, 0.8383293268314073, 0.03158184731104894, 0.005079517400465046, 0.09807639052204041, 0.00021267104441893782, 0.002812435050440354, 0.09905452804226003]
 ⋮
 [0.07298745597949036, 0.24897247623820626, 0.8383007742866869, 0.03158103041201427, 0.005079478820963999, 0.09807399035677364, 0.0002127329960413663, 0.0028123966020952708, 0.09905395275351637]
 [0.0729874778352529, 0.2489724877765076, 0.8382973953218201, 0.031580831980856636, 0.005079474232119774, 0.09807372554833006, 0.00021274107590946308, 0.0028123912202478547, 0.09905384529873094]
 [0.07298750043744294, 0.2489724953320386, 0.8382941760341811, 0.031580633648816324, 0.005079469861346299, 0.09807347528175157, 0.0002127488654749285, 0.002812385997122885, 0.09905373820373098]
 [0.07298752361226937, 0.2489724983136289, 0.8382911510573068, 0.031580437058832296, 0.005079465756123402, 0.09807324256173498, 0.00021275628614032864, 0.002812380983238756, 0.09905363237897294]
 [0.07298754718594118, 0.24897249613010797, 0.8382883550247344, 0.03158024385384352, 0.005079461963930914, 0.09807303039297703, 0.0002127632593082294, 0.002812376229113863, 0.09905352873491317]
 [0.07298757098466731, 0.2489724881903054, 0.8382858225700007, 0.03158005567678893, 0.005079458532248665, 0.09807284178017452, 0.0002127697063811969, 0.0028123717852666006, 0.09905342818200816]
 [0.0729875948346567, 0.2489724739030508, 0.8382835883266426, 0.03157987417060752, 0.005079455508556486, 0.09807267972802418, 0.00021277554876179713, 0.0028123677022153628, 0.09905333163071424]
 [0.07298761856211827, 0.2489724526771737, 0.8382816869281972, 0.031579700978238215, 0.005079452940334206, 0.09807254724122283, 0.00021278070785259612, 0.0028123640304785454, 0.09905323999148788]
 [0.07298764199326101, 0.24897242392150368, 0.8382801530082012, 0.03157953774262, 0.005079450875061657, 0.09807244732446721, 0.00021278510505615987, 0.0028123608205745413, 0.09905315417478545]

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/fab42f6178325bfc548ac92fe66c1517a406c85a9b9774492a304acf59f439c1.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, [Glc => 10mM, kATPCa=>10Hz/mM, kATP=>0.055Hz])
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.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/hlGut/src/deprecations.jl:45
_images/e2daf25ff010417fc7de719c8b394593d081303c4bcf2afe5212286db308c996.png

kATPCa : 90 -> 0.1

prob4 = ODEProblem(sysosci, [], tend, [Glc => 10mM, kATPCa=>0.1Hz/mM, kATP=>0.06Hz])
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.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/hlGut/src/deprecations.jl:45
_images/02f6363acca36af9fb7490add09bf6abbb431e1b62b97dd3880c644b1fb806f2.png

This notebook was generated using Literate.jl.