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}:
      2.196435053752864
      0.0007273589705197244
      0.000981691808598968
    245.22202099807583
    281.0719264404058
    197.06621753790645
     34.41477140951692
     60.37913837798986
    229.2919970960533
     22.321405235243013
      ⋮
     23.148552252520883
      0.6928662120941113
      0.20648605639609985
   1295.5364065086924
   1299.7940884726063
      0.20381438659933854
 145551.11317854663
  10276.480036234318
    -85.29107641956527

@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

 20.469907 seconds (8.28 M allocations: 3.680 GiB, 2.61% gc time, 12.24% compilation time)
sys.gssg_i => 0.6266485048629427,
sys.h2o2_i => 0.0002160109305120311,
sys.sox_i => 0.00028193442352659245,
sys.cytc_ox => 190.99032016698382,
sys.cytc1_ox => 243.39273317769164,
sys.fes_ox => 138.12099168741568,
sys.blo_bhr => 25.693323786169316,
sys.blr_bho => 67.38311339022796,
sys.blo_bho => 231.22615960469963,
sys.QH2_p => 0.5304090644069537,
sys.QH2_n => 0.5681269384558159,
sys.Q_n => 1972.1657466314166,
sys.SQn => 54.11028700204462,
sys.SQp => 0.4219759224046464,
sys.N2r_C1 => 16.962752913631025,
sys.oaa => 0.0034565339290617797,
sys.mal => 17.461777311088888,
sys.fum => 30.284280329284,
sys.suc => 0.7824266335897897,
sys.scoa => 70.28220061404784,
sys.akg => 0.4837564931870146,
sys.isoc => 808.1303501210984,
sys.j_na => 0.9785813135724989,
sys.h_na => 0.9700523872313703,
sys.m_na => 0.001453690150356163,
sys.x_k => 0.0028256803347640196,
sys.htr_ca => 137.6817784377457,
sys.ltr_ca => 23.565903213334103,
sys.x_n1 => 0.020515032933281677,
sys.x_p3 => 0.010594532589095668,
sys.x_p2 => 0.01214123098542219,
sys.x_p1 => 0.006469843750282456,
sys.x_p0 => 0.00727590633038226,
sys.crp_ic => 16677.260291033977,
sys.crp_i => 16747.116810150023,
sys.adp_ic => 37.74866681710288,
sys.x_yca => 0.7245412129389334,
sys.cca4_lcc => 9.806251992166545e-23,
sys.cca3_lcc => 3.119442863967385e-17,
sys.cca2_lcc => 3.721178309622979e-12,
sys.cca1_lcc => 1.9728745364249864e-7,
sys.cca0_lcc => 0.00392235910123243,
sys.o_lcc => 1.4711284914459437e-23,
sys.c4_lcc => 9.727492988731676e-23,
sys.c3_lcc => 1.2377050359007134e-16,
sys.c2_lcc => 5.90588827911184e-11,
sys.c1_lcc => 1.2524799881224833e-5,
sys.pc2_ryr => 0.12682911264081026,
sys.po2_ryr => 2.8084081198505617e-9,
sys.po1_ryr => 0.0001725424697079885,
sys.O2 => 5.920152808057068,
sys.sox_m => 0.11424829460620131,
sys.dpsi => 171.47543371889324,
sys.nadh_m => 441.0285389318614,
sys.adp_m => 6.224895463254149,
sys.adp_i => 37.12036497085859,
sys.ca_m => 0.851767834462317,
sys.ca_ss => 0.1979210800388898,
sys.ca_jsr => 1238.368871750374,
sys.ca_nsr => 1242.3601900850263,
sys.ca_i => 0.19518420240776727,
sys.k_i => 144972.33082726505,
sys.na_i => 10213.483972685903,
sys.vm => -85.36727678461618,

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.