Initial conditions

1Hz pacing for 1000 seconds.

using OrdinaryDiffEq
using ModelingToolkit
using ECMEDox
using ECMEDox: second, mM, Hz, μM, mV
using Plots
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; callback = callback, saveat=0.01second);
 21.841534 seconds (20.89 M allocations: 990.471 MiB, 0.61% gc time, 36.19% compilation time)
for i in sts
    istr = replace(string(i), "(t)" => "")
    println("sys.", istr, " => ", sol[i][end], ",")
end
sys.gssg_i => 0.6265895252239285,
sys.h2o2_i => 0.00021599097781387641,
sys.sox_i => 0.000281908019349653,
sys.cytc_ox => 190.9942064650482,
sys.cytc1_ox => 243.39572775560168,
sys.fes_ox => 138.1247054108801,
sys.blo_bhr => 25.693443738402454,
sys.blr_bho => 67.38162764992859,
sys.blo_bho => 231.22754320816577,
sys.QH2_p => 0.5303653651617399,
sys.QH2_n => 0.5680839234749541,
sys.Q_n => 1972.1657215559703,
sys.SQn => 54.11043560924573,
sys.SQp => 0.4219634954914862,
sys.N2r_C1 => 16.96274764971865,
sys.oaa => 0.0034565870103478637,
sys.mal => 17.459086815466467,
sys.fum => 30.281799654859434,
sys.suc => 0.782440772356075,
sys.scoa => 70.27336843756892,
sys.akg => 0.48375017081405797,
sys.isoc => 808.1399060773019,
sys.j_na => 0.9785809587363049,
sys.h_na => 0.9700517570637847,
sys.m_na => 0.0014537117796979628,
sys.x_k => 0.0028261999315035095,
sys.htr_ca => 137.68179458376923,
sys.ltr_ca => 23.566070750638414,
sys.x_n1 => 0.020516029261279715,
sys.x_p3 => 0.010595307782823218,
sys.x_p2 => 0.012142111983808474,
sys.x_p1 => 0.006470302862256833,
sys.x_p0 => 0.00727635236662959,
sys.crp_ic => 16676.91890175737,
sys.crp_i => 16746.777783513873,
sys.adp_ic => 37.75097704512865,
sys.x_yca => 0.7245465013652808,
sys.cca4_lcc => 9.809958217150686e-23,
sys.cca3_lcc => 3.120039381514477e-17,
sys.cca2_lcc => 3.721593607602344e-12,
sys.cca1_lcc => 1.9730128695085108e-7,
sys.cca0_lcc => 0.003922560427831106,
sys.o_lcc => 2.0665952349622882e-23,
sys.c4_lcc => 9.739899479051778e-23,
sys.c3_lcc => 1.2377646326795985e-16,
sys.c2_lcc => 5.90607417504726e-11,
sys.c1_lcc => 1.2524994845841358e-5,
sys.pc2_ryr => 0.1268315869668816,
sys.po2_ryr => 2.8085358346184804e-9,
sys.po1_ryr => 0.0001725460552723041,
sys.O2 => 5.920151510561289,
sys.sox_m => 0.11424499555126147,
sys.dpsi => 171.4740053680407,
sys.nadh_m => 440.91731505543635,
sys.adp_m => 6.2257878819562915,
sys.adp_i => 37.12263459692962,
sys.ca_m => 0.8518010873547079,
sys.ca_ss => 0.19792262061120708,
sys.ca_jsr => 1238.3749208769589,
sys.ca_nsr => 1242.3665362878057,
sys.ca_i => 0.1951856729070031,
sys.k_i => 144972.08735520177,
sys.na_i => 10213.437817329735,
sys.vm => -85.36718765864508,

Action potential

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

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")

CAC flux

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

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)")

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

Proton pumping

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

ROS

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

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

plot(sol, idxs=100 * sys.vROS / (sys.vO2 + sys.vROS), title="O2 Shunt", tspan=(900e3, 910e3), ylims=(0, 5))

MMP

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

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

Runtime information

using InteractiveUtils
InteractiveUtils.versioninfo()
Julia Version 1.12.4
Commit 01a2eadb047 (2026-01-06 16:56 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 4 × AMD EPYC 7763 64-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, znver3)
  GC: Built with stock GC
Threads: 4 default, 1 interactive, 4 GC (on 4 virtual cores)
Environment:
  JULIA_CPU_TARGET = generic;icelake-server,clone_all;znver3,clone_all
  JULIA_CONDAPKG_OFFLINE = true
  JULIA_CONDAPKG_BACKEND = Null
  JULIA_CI = true
  LD_LIBRARY_PATH = /opt/hostedtoolcache/Python/3.14.6/x64/lib
  JULIA_NUM_THREADS = auto
using Pkg
Pkg.status()
Project ECMEDox v0.4.0
Status `~/work/ecme-rirr-dox/ecme-rirr-dox/Project.toml`
⌃ [459566f4] DiffEqCallbacks v4.11.0
  [0b91fe84] DisplayAs v0.1.6
⌃ [f6369f11] ForwardDiff v1.3.1
⌃ [682c06a0] JSON v1.3.0
⌃ [23fbe1c1] Latexify v0.16.10
  [98b081ad] Literate v2.21.0
⌃ [961ee093] ModelingToolkit v10.31.2
⌃ [77ba4419] NaNMath v1.1.3
⌃ [8913a72c] NonlinearSolve v4.13.0
⌃ [1dea7af3] OrdinaryDiffEq v6.105.0
⌃ [91a5bcdd] Plots v1.41.4
  [33c8b6b6] ProgressLogging v0.1.6
⌃ [9672c7b4] SteadyStateDiffEq v2.9.0
⌅ [0c5d862f] Symbolics v6.58.0
  [37e2e46d] LinearAlgebra v1.12.0
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`

This notebook was generated using Literate.jl.

Back to top