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

 33.774092 seconds (17.55 M allocations: 969.823 MiB, 1.50% gc time, 25.72% compilation time)
sys.gssg_i => 0.6266485048275071,
sys.h2o2_i => 0.00021601093049525123,
sys.sox_i => 0.00028193442350450584,
sys.cytc_ox => 190.9903201686805,
sys.cytc1_ox => 243.39273317895325,
sys.fes_ox => 138.1209916884662,
sys.blo_bhr => 25.69332378721287,
sys.blr_bho => 67.38311339081238,
sys.blo_bho => 231.22615960303426,
sys.QH2_p => 0.530409064403493,
sys.QH2_n => 0.568126938455307,
sys.Q_n => 1972.1657466307868,
sys.SQn => 54.11028700330599,
sys.SQp => 0.4219759224039639,
sys.N2r_C1 => 16.962752913626485,
sys.oaa => 0.0034565339292543375,
sys.mal => 17.46177730957997,
sys.fum => 30.28428032846593,
sys.suc => 0.7824266336379635,
sys.scoa => 70.28220060253133,
sys.akg => 0.483756492978476,
sys.isoc => 808.1303501304549,
sys.j_na => 0.9785813134557223,
sys.h_na => 0.9700523870235357,
sys.m_na => 0.0014536901576709923,
sys.x_k => 0.002825680345259897,
sys.htr_ca => 137.68177843818987,
sys.ltr_ca => 23.565903216891847,
sys.x_n1 => 0.020515032950018264,
sys.x_p3 => 0.010594532597918219,
sys.x_p2 => 0.012141230995523885,
sys.x_p1 => 0.0064698437556523014,
sys.x_p0 => 0.00727590633630333,
sys.crp_ic => 16677.260290265833,
sys.crp_i => 16747.116809365325,
sys.adp_ic => 37.748666822300905,
sys.x_yca => 0.7245412128492645,
sys.cca4_lcc => 9.806252163775228e-23,
sys.cca3_lcc => 3.119442914709192e-17,
sys.cca2_lcc => 3.72117835025185e-12,
sys.cca1_lcc => 1.9728745474243374e-7,
sys.cca0_lcc => 0.003922359102155739,
sys.o_lcc => 1.4411487478248164e-23,
sys.c4_lcc => 9.726953877075017e-23,
sys.c3_lcc => 1.2377050556458878e-16,
sys.c2_lcc => 5.905888342168872e-11,
sys.c1_lcc => 1.252479994808198e-5,
sys.pc2_ryr => 0.12682911261082286,
sys.po2_ryr => 2.80840812120647e-9,
sys.po1_ryr => 0.00017254246969060553,
sys.O2 => 5.92015280805098,
sys.sox_m => 0.11424829460595995,
sys.dpsi => 171.47543371718447,
sys.nadh_m => 441.02853883784786,
sys.adp_m => 6.224895464638173,
sys.adp_i => 37.120364976127334,
sys.ca_m => 0.8517678350226103,
sys.ca_ss => 0.1979210800774621,
sys.ca_jsr => 1238.3688720650182,
sys.ca_nsr => 1242.3601904012212,
sys.ca_i => 0.19518420244591997,
sys.k_i => 144972.3306562517,
sys.na_i => 10213.483968167471,
sys.vm => -85.36727675443937,

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.