Initial conditions

1Hz pacing for 1000 seconds

using ProgressLogging
using OrdinaryDiffEq
using ModelingToolkit
using ECMEDox
using ECMEDox: second, mM, Hz, μM, mV
using Plots
using DisplayAs: PNG

tend = 1000.0second
bcl = 1.0second
@named sys = build_model()
u0 = build_u0(sys)
sts = unknowns(sys)
alg = KenCarp47()
@unpack iStim = sys
callback = build_stim_callbacks(iStim, tend; period=bcl)
prob = ODEProblem(sys, u0, tend)

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Initialization status: FULLY_DETERMINED
Non-trivial mass matrix: false
timespan: (0.0, 1.0e6)
u0: 64-element Vector{Float64}:
      1.4355864646916074
      0.0004949419698468167
      0.0006578256307800626
    209.31680906407098
    257.044146799667
    156.82453780216784
     54.88374367717081
    128.68082835522182
    136.13456850129444
    168.9705244901982
      ⋮
     43.96538500887081
      0.8926956182412644
      0.19547857075761974
   1217.9822251524322
   1222.0613051409616
      0.19267665443392065
 144821.2542377956
  10200.108762034
    -85.39304888131878

@time sol = solve(prob, alg; reltol=1e-6, abstol=1e-6, progress=true, callback = callback)

for i in sts
    istr = replace(string(i), "(t)" => "")
    println("sys.", istr, " => ", sol[i][end], ",")
end

 16.897973 seconds (17.62 M allocations: 1013.964 MiB, 1.71% gc time, 30.82% compilation time)
sys.gssg_i => 1.4935057071642388,
sys.h2o2_i => 0.0005136979121467542,
sys.sox_i => 0.0006835947826241502,
sys.cytc_ox => 208.84377105075984,
sys.cytc1_ox => 256.7039416971413,
sys.fes_ox => 156.32093110708306,
sys.blo_bhr => 55.232756210612884,
sys.blr_bho => 129.87267376143038,
sys.blo_bho => 134.44204128949877,
sys.QH2_p => 174.51398146668117,
sys.QH2_n => 174.55189407761645,
sys.Q_n => 1748.7946569305636,
sys.SQn => 152.27842137727842,
sys.SQp => 1.0284813447942511,
sys.N2r_C1 => 16.964060683838348,
sys.oaa => 0.0034681157233715835,
sys.mal => 18.29946283716348,
sys.fum => 31.180206035228736,
sys.suc => 0.8585527996607384,
sys.scoa => 40.601230416373696,
sys.akg => 0.4689786821865219,
sys.isoc => 827.3397731879196,
sys.j_na => 0.9786768078012272,
sys.h_na => 0.9702225491707038,
sys.m_na => 0.0014476950418653598,
sys.x_k => 0.0028087252421126,
sys.x_n1 => 0.01969120424061365,
sys.x_p3 => 0.010100638687329737,
sys.x_p2 => 0.011574355321828777,
sys.x_p1 => 0.006166457736490283,
sys.x_p0 => 0.0069322729407981214,
sys.htr_ca => 137.6535306904159,
sys.ltr_ca => 23.345201026541336,
sys.crp_ic => 15743.286269805632,
sys.crp_i => 15814.437735125637,
sys.adp_ic => 44.43819163884076,
sys.pc2_ryr => 0.13358010445440546,
sys.po2_ryr => 2.8153712809619514e-9,
sys.po1_ryr => 0.00017932762230896102,
sys.x_yca => 0.724768808214493,
sys.cca4_lcc => 9.53208479923466e-23,
sys.cca3_lcc => 3.045554944440158e-17,
sys.cca2_lcc => 3.64900325321118e-12,
sys.cca1_lcc => 1.9431102711312328e-7,
sys.cca0_lcc => 0.00388015812388709,
sys.o_lcc => 1.4327054738738193e-23,
sys.c4_lcc => 9.558344530823941e-23,
sys.c3_lcc => 1.2215828884361993e-16,
sys.c2_lcc => 5.854573000754625e-11,
sys.c1_lcc => 1.2470532753689835e-5,
sys.O2 => 5.918838279399363,
sys.sox_m => 0.2774985950519481,
sys.dpsi => 165.66086399182106,
sys.nadh_m => 469.8008506298959,
sys.adp_m => 10.847572882371194,
sys.adp_i => 43.731105943801765,
sys.ca_m => 0.8917084518081624,
sys.ca_ss => 0.19555352504252876,
sys.ca_jsr => 1218.6732697138798,
sys.ca_nsr => 1222.7471749179485,
sys.ca_i => 0.19275426470601623,
sys.k_i => 144824.80761399592,
sys.na_i => 10200.467630348216,
sys.vm => -85.39205943050112,

plot(sol, idxs=sys.vm, legend=:right, tspan=(900second, 901second)) |> PNG

Citric acid cycle metabolites Citrate and isocitrate have the highest concentrations

@unpack cit, isoc, oaa, akg, scoa, suc, fum, mal = sys
plot(sol, idxs=[cit, isoc, oaa, akg, scoa, suc, fum, mal], legend=:right, title="CAC metabolites") |> PNG

CAC flux

@unpack vMDH, vAAT, vIDH = sys
plot(sol, idxs=[vMDH, vAAT, vIDH], legend=:right, title="CAC flux") |> PNG

Q cycle

@unpack Q_n, SQn, QH2_n, QH2_p, Q_p, SQp, fes_ox, fes_rd, cytc_ox, cytc_rd = sys
plot(sol, idxs=[Q_n + Q_p, SQn, QH2_n + QH2_p, SQp], title="Q cycle", legend=:left, xlabel="Time (ms)", ylabel="Conc. (μM)") |> PNG

Q cycle downstream

plot(sol, idxs=[fes_ox, fes_rd, cytc_ox, cytc_rd], title="Q cycle (downstream)", legend=:left, xlabel="Time (ms)", ylabel="Conc. (μM)") |> PNG

plot(sol, idxs = [sys.vHresC1, sys.vHresC3, sys.vHresC4], ylims=(0, 3)) |> PNG

plot(sol, idxs = [sys.sox_i, sys.sox_m], tspan=(900e3, 910e3)) |> PNG

ROS

plot(sol, idxs = [sys.vROSIf, sys.vROSIq, sys.vROSC1, sys.vROSC3], tspan=(900e3, 910e3)) |> PNG

O2 Shunt

plot(sol, idxs=100 * sys.vROS / (sys.vO2 + sys.vROS), title="O2 Shunt", tspan=(900e3, 910e3)) |> PNG

MMP

plot(sol, idxs = [sys.dpsi], tspan=(900e3, 910e3)) |> PNG

ATP synthesis rate

plot(sol, idxs = [sys.vC5], tspan=(900e3, 910e3)) |> PNG

---

This notebook was generated using Literate.jl.