using OrdinaryDiffEq
using ModelingToolkit
using MitochondrialDynamics
using MitochondrialDynamics: second, μM, mV, mM, Hz, minute, hil
import PythonPlot as plt
plt.matplotlib.rcParams["font.size"] = 1414@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)24-element Vector{Symbolics.Equation}:
x1(t) ~ 1 - 3x3(t) - 2x2(t)
ATP_c(t) ~ ((1 - NaNMath.sqrt(1 + 4(-1 + 4kEqAK)*(1 - AEC(t))*AEC(t))) / (-2 + 8kEqAK) + AEC(t))*ΣAc
ADP_c(t) ~ (2(-1 + NaNMath.sqrt(1 + 4(-1 + 4kEqAK)*(1 - AEC(t))*AEC(t)))*ΣAc) / (-2 + 8kEqAK)
J_HL(t) ~ pHleak*rHL*exp(kvHleak*ΔΨm(t))
NAD_c(t) ~ -NADH_c(t) + Σn_c
NAD_m(t) ~ -NADH_m(t) + Σn_m
J_HR(t) ~ ((1 + KaETC*ΔΨm(t))*VmaxETC*rETC*NADH_m(t)) / ((1 + KbETC*ΔΨm(t))*(KnadhETC + NADH_m(t)))
degavg(t) ~ (3x3(t) + 2x2(t) + x1(t)) / (x3(t) + x2(t) + x1(t))
J_GK(t) ~ (VmaxGK*ATP_c(t)*NaNMath.pow(Glc, nGK)) / ((KatpGK + ATP_c(t))*(NaNMath.pow(KglcGK, nGK) + NaNMath.pow(Glc, nGK)))
Ca_c(t) ~ RestingCa - ((-ActivatedCa*(1 + 0.8(-1 + 1250(exp(1 - 5rem(t / 120.0, 1.0))^4)*(rem(t / 120.0, 1.0)^4)))*NaNMath.pow(ATP_c(t), NCac)) / (NaNMath.pow(ATP_c(t), NCac) + NaNMath.pow(KatpCac*ADP_c(t), NCac)))
⋮
J_PDH(t) ~ (VmaxPDH*rPDH*NAD_m(t)*Pyr(t)) / ((NAD_m(t) + KnadPDH*NADH_m(t))*(KpyrPDH + Pyr(t))*(1 + (1 + U1PDH*((KcaPDH / (KcaPDH + Ca_m(t)))^2))*U2PDH))
J_NADHT(t) ~ (VmaxNADHT*NADH_c(t)*NAD_m(t)) / ((NADH_c(t) + Ktn_c*NAD_c(t))*(NAD_m(t) + Ktn_m*NADH_m(t)))
J_O2(t) ~ (1//10)*J_HR(t)
J_NCLX(t) ~ (VmaxNCLX*((-((Na_m / KnaNCLX)^2)*Ca_c(t)) / KcaNCLX + (Ca_m(t)*((Na_c / KnaNCLX)^2)) / KcaNCLX)) / (1 + (((Na_m / KnaNCLX)^2)*Ca_c(t)) / KcaNCLX + Ca_m(t) / KcaNCLX + (Ca_m(t)*((Na_c / KnaNCLX)^2)) / KcaNCLX + Ca_c(t) / KcaNCLX + (Na_m / KnaNCLX)^2 + (Na_c / KnaNCLX)^2)
J_MCU(t) ~ (74.87176701560523PcaMCU*(-0.2Ca_m(t) + 0.341Ca_c(t)*(1 + expm1(74.87176701560523ΔΨm(t))))*ΔΨm(t)) / expm1(74.87176701560523ΔΨm(t))
J_ANT(t) ~ (1//3)*J_HF(t)
J_DH(t) ~ J_FFA(t) + 4.6J_PDH(t)
tiptip(t) ~ -((-kfuse1*J_ANT(t)*(x1(t)^2)) / J_HL(t)) - kfiss1*x2(t)
tipside(t) ~ -((-kfuse2*J_ANT(t)*x2(t)*x1(t)) / J_HL(t)) - kfiss2*x3(t)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.
└ @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/nrHl8/src/deprecations.jl:53
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.030896192894865025, 0.21295077836820114, 0.9825638030195014, 0.03474219851893958, 0.005308721776598269, 0.12217147941950664, 0.0008340349530882978, 0.0031096559823532415, 0.10663891083237557]
[0.030982379476767317, 0.21323952995228992, 0.9835907617419575, 0.03492136668151196, 0.0053139165415342715, 0.12422676423355648, 0.0008280947235631401, 0.0031206800661563636, 0.10714795524026766]
[0.031046708124656815, 0.21344994205959628, 0.9845817541980271, 0.03514386763699833, 0.005319553127195548, 0.12654195778428817, 0.0008189357919522203, 0.003133866105919805, 0.10776555995464054]
[0.031089305175396164, 0.21358312541178587, 0.9855064562638932, 0.03540844553019241, 0.005325439574049158, 0.12904032393241907, 0.0008065021071712109, 0.0031489859856371314, 0.10848050649979553]
[0.031110621670136094, 0.2136410848116792, 0.9863455302510558, 0.03571324296841961, 0.005331398493786738, 0.13164496724824126, 0.0007907654204095447, 0.003165765260230514, 0.10927892417712361]
[0.031112395812853424, 0.2136305076195222, 0.9870623310628815, 0.036050020869946554, 0.005337034660269277, 0.13417595013308653, 0.0007719416439220301, 0.0031833889152117403, 0.11012063689368846]
[0.031098719578800296, 0.21356635181183706, 0.9876189932805626, 0.03639819460775611, 0.005341781237930851, 0.13636356762887938, 0.0007507952400103626, 0.0032002697237616846, 0.11092935371456224]
[0.03107135984497027, 0.21345484990143274, 0.9880367260510318, 0.0367480897427108, 0.005345672120024203, 0.13821441993653236, 0.0007277429081124214, 0.0032159946206352396, 0.11168725819271969]
[0.031032745029412164, 0.213304527236694, 0.988348134701614, 0.03708683521861731, 0.005348797780151491, 0.13975221961789466, 0.0007034742232345601, 0.0032301262775251017, 0.11237423519655106]
[0.030986769188265678, 0.21312901409343388, 0.9885970389628247, 0.03739637851060485, 0.0053513199539413185, 0.14101771078569325, 0.000679341909715767, 0.003242278942731089, 0.11296874668406232]
⋮
[0.029699206579523997, 0.20863909504947872, 0.9791205090970783, 0.03503452437907707, 0.0052945060453428285, 0.11691636512073889, 0.0007580585505795374, 0.0030932605002094416, 0.10601958149508883]
[0.029868087120471596, 0.20927347297678617, 0.978715177024182, 0.03472996293439517, 0.005292834883146665, 0.11629476217235214, 0.0007841068095212187, 0.003086433731169731, 0.10567378053197109]
[0.030040318127384693, 0.20991165179824006, 0.9786647804746441, 0.03454825665506743, 0.005292536520063341, 0.11617222655113119, 0.0008035761608557764, 0.003083293877430465, 0.10550395631286293]
[0.03021019315541264, 0.21053282675706458, 0.9788889260326835, 0.03445403805393275, 0.005293291425341006, 0.11643122771230659, 0.0008177967249617406, 0.0030828885877388005, 0.10546527576492619]
[0.03037308976050629, 0.2111206796406811, 0.9793252974621159, 0.034421162149261886, 0.005294861694270066, 0.11698109696720943, 0.0008278184692093362, 0.00308453612385961, 0.10552286032834006]
[0.030526676790023206, 0.21166750996040276, 0.9799548119794207, 0.034441181537933775, 0.005297244418379861, 0.1178336709784703, 0.0008340182270240145, 0.003088101614582709, 0.10567196519389715]
[0.03066608755464715, 0.2121578180551516, 0.980718753237448, 0.03450208288266646, 0.005300360263680669, 0.11897989040071114, 0.0008370682453207643, 0.0030934500386995637, 0.10590595222622347]
[0.03078949393860513, 0.2125860410469203, 0.9815961903119094, 0.034602160601498975, 0.005304190203655462, 0.12042265165529296, 0.0008370965073495683, 0.0031006031389581223, 0.10622633504146481]
[0.030895600444306233, 0.21294832240577266, 0.9825733983421877, 0.03474115375904537, 0.005308727660128024, 0.1221677321926966, 0.0008341470122119232, 0.0031096007391064138, 0.10663570308182106]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
endplot_fig5 (generic function with 2 methods)fig5 = plot_fig5(sol)
Export figure
exportTIF(fig5, "Fig6-ca-oscillation.tif")Python: NoneTuning 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.
└ @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/nrHl8/src/deprecations.jl:53
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.
└ @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/nrHl8/src/deprecations.jl:53
This notebook was generated using Literate.jl.