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

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}, 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!, 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{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{NAD\_m}\left( t \right) &= - \mathtt{NADH\_m}\left( t \right) + \mathtt{{\Sigma}n\_m} \\ \mathtt{NAD\_c}\left( t \right) &= - \mathtt{NADH\_c}\left( t \right) + \mathtt{{\Sigma}n\_c} \\ \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{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{J\_HL}\left( t \right) &= \mathtt{pHleak} \mathtt{rHL} e^{\mathtt{kvHleak} \mathtt{\Delta{\Psi}m}\left( t \right)} \\ \mathtt{x1}\left( t \right) &= 1 - 3 \mathtt{x3}\left( t \right) - 2 \mathtt{x2}\left( t \right) \\ \mathtt{J\_O2}\left( t \right) &= \frac{1}{10} \mathtt{J\_HR}\left( t \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\_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\_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{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\_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{AMP\_c}\left( t \right) &= - \mathtt{ATP\_c}\left( t \right) - \mathtt{ADP\_c}\left( t \right) + \mathtt{{\Sigma}Ac} \\ \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\_DH}\left( t \right) &= \mathtt{J\_FFA}\left( t \right) + 4.6 \mathtt{J\_PDH}\left( t \right) \\ \mathtt{J\_ANT}\left( t \right) &= \frac{1}{3} \mathtt{J\_HF}\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{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) \\ \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) \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: The initialization system is structurally singular. Guess values may significantly affect the initial values of the ODE. The problematic variables are Any[J_MCU(t)].
│ 
│ Note that the identification of problematic variables is a best-effort heuristic.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/vIsAP/src/systems/diffeqs/abstractodesystem.jl:1426
┌ Warning: Initialization system is underdetermined. 2 equations for 3 unknowns. Initialization will default to using least squares. `SCCNonlinearProblem` can only be used for initialization of fully determined systems and hence will not be used here. To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/vIsAP/src/systems/diffeqs/abstractodesystem.jl:1456

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.09905459455316733, 0.0028124343488794986, 0.0002126738642724931, 0.0980763352877003, 0.005079515343908573, 0.03158188299843019, 0.8383275916263324, 0.24897140514600974, 0.07298768624306845]
 [0.0990545766220652, 0.00281243378404554, 0.00021267452146499937, 0.09807630785119426, 0.005079514977486804, 0.03158186074028628, 0.8383273436364552, 0.2489713650818988, 0.07298765519282241]
 [0.09905456148956045, 0.002812433397013501, 0.00021267490775616924, 0.09807628978868912, 0.005079514778174247, 0.031581843441071594, 0.8383272196614613, 0.24897132715525125, 0.07298762486822352]
 [0.0990545490856472, 0.0028124331804523987, 0.00021267503526340738, 0.09807628067993324, 0.0050795147380949535, 0.03158183093548211, 0.8383272137796621, 0.2489712913825985, 0.07298759525189745]
 [0.0990545393403196, 0.0028124331270312493, 0.0002126749161041186, 0.09807628010467498, 0.005079514849372985, 0.0315818230582138, 0.8383273200693696, 0.2489712577804719, 0.07298756632646983]
 [0.09905453218357173, 0.0028124332294190675, 0.00021267456239570763, 0.09807628764266267, 0.0050795151041324, 0.031581819643962616, 0.8383275326088948, 0.2489712263654028, 0.07298753807456629]
 [0.09905452754539772, 0.0028124334802848714, 0.0002126739862555793, 0.09807630287364472, 0.0050795154944972545, 0.031581820527424566, 0.8383278454765496, 0.24897119715392263, 0.07298751047881245]
 [0.09905452535579169, 0.002812433872297676, 0.0002126731998011384, 0.09807632537736945, 0.005079516012591609, 0.03158182554329562, 0.8383282527506456, 0.24897117016256273, 0.072987483521834]
 [0.09905452554474775, 0.002812434398126498, 0.00021267221514978963, 0.09807635473358522, 0.00507951665053952, 0.03158183452627175, 0.8383287485094941, 0.24897114540785456, 0.07298745718625652]
 [0.09905452804226005, 0.0028124350504403536, 0.00021267104441893776, 0.09807639052204042, 0.005079517400465046, 0.031581847311048926, 0.8383293268314069, 0.24897112290632942, 0.07298743145470567]
 ⋮
 [0.09905395275351637, 0.002812396602095271, 0.00021273299604136628, 0.09807399035677364, 0.005079478820963998, 0.03158103041201425, 0.8383007742866869, 0.24897247623820626, 0.07298745597949037]
 [0.09905384529873094, 0.0028123912202478547, 0.00021274107590946308, 0.09807372554833005, 0.005079474232119774, 0.031580831980856615, 0.83829739532182, 0.24897248777650763, 0.0729874778352529]
 [0.09905373820373099, 0.0028123859971228853, 0.0002127488654749285, 0.09807347528175157, 0.005079469861346298, 0.0315806336488163, 0.8382941760341811, 0.24897249533203866, 0.07298750043744294]
 [0.09905363237897294, 0.0028123809832387566, 0.00021275628614032862, 0.09807324256173497, 0.005079465756123401, 0.03158043705883228, 0.8382911510573068, 0.24897249831362891, 0.07298752361226939]
 [0.0990535287349132, 0.002812376229113864, 0.00021276325930822938, 0.09807303039297702, 0.005079461963930914, 0.031580243853843504, 0.8382883550247344, 0.24897249613010797, 0.07298754718594119]
 [0.09905342818200816, 0.0028123717852666015, 0.00021276970638119688, 0.0980728417801745, 0.005079458532248665, 0.03158005567678892, 0.8382858225700007, 0.24897248819030543, 0.07298757098466732]
 [0.09905333163071427, 0.002812367702215364, 0.00021277554876179713, 0.09807267972802418, 0.005079455508556486, 0.0315798741706075, 0.8382835883266426, 0.24897247390305083, 0.0729875948346567]
 [0.09905323999148789, 0.0028123640304785463, 0.0002127807078525961, 0.09807254724122283, 0.005079452940334206, 0.0315797009782382, 0.8382816869281972, 0.24897245267717372, 0.07298761856211829]
 [0.09905315417478548, 0.0028123608205745426, 0.00021278510505615985, 0.09807244732446721, 0.005079450875061657, 0.031579537742619984, 0.8382801530082012, 0.2489724239215037, 0.07298764199326102]

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: The initialization system is structurally singular. Guess values may significantly affect the initial values of the ODE. The problematic variables are Any[J_MCU(t)].
│ 
│ Note that the identification of problematic variables is a best-effort heuristic.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/vIsAP/src/systems/diffeqs/abstractodesystem.jl:1426
┌ Warning: Initialization system is underdetermined. 2 equations for 3 unknowns. Initialization will default to using least squares. `SCCNonlinearProblem` can only be used for initialization of fully determined systems and hence will not be used here. To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/vIsAP/src/systems/diffeqs/abstractodesystem.jl:1456
_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: The initialization system is structurally singular. Guess values may significantly affect the initial values of the ODE. The problematic variables are Any[J_MCU(t)].
│ 
│ Note that the identification of problematic variables is a best-effort heuristic.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/vIsAP/src/systems/diffeqs/abstractodesystem.jl:1426
┌ Warning: Initialization system is underdetermined. 2 equations for 3 unknowns. Initialization will default to using least squares. `SCCNonlinearProblem` can only be used for initialization of fully determined systems and hence will not be used here. To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/vIsAP/src/systems/diffeqs/abstractodesystem.jl:1456
_images/02f6363acca36af9fb7490add09bf6abbb431e1b62b97dd3880c644b1fb806f2.png

This notebook was generated using Literate.jl.