Initial conditions

1Hz pacing for 1000 seconds.

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; callback = callback, saveat=0.01second);
 21.812590 seconds (20.85 M allocations: 988.208 MiB, 1.42% gc time, 39.51% compilation time)
for i in sts
    istr = replace(string(i), "(t)" => "")
    println("sys.", istr, " => ", sol[i][end], ",")
end
sys.gssg_i => 0.6265895252227057,
sys.h2o2_i => 0.0002159909778123485,
sys.sox_i => 0.00028190801934770745,
sys.cytc_ox => 190.9942064651759,
sys.cytc1_ox => 243.39572775569582,
sys.fes_ox => 138.1247054109357,
sys.blo_bhr => 25.69344373847062,
sys.blr_bho => 67.38162764998943,
sys.blo_bho => 231.22754320803372,
sys.QH2_p => 0.5303653651619983,
sys.QH2_n => 0.5680839234754813,
sys.Q_n => 1972.1657215559387,
sys.SQn => 54.11043560930814,
sys.SQp => 0.4219634954912853,
sys.N2r_C1 => 16.962747649718544,
sys.oaa => 0.00345658701035074,
sys.mal => 17.459086815402173,
sys.fum => 30.281799654823033,
sys.suc => 0.7824407723640625,
sys.scoa => 70.27336843732768,
sys.akg => 0.48375017081403227,
sys.isoc => 808.1399060775319,
sys.j_na => 0.978580958736307,
sys.h_na => 0.9700517570637879,
sys.m_na => 0.0014537117796978982,
sys.x_k => 0.0028261999315025966,
sys.htr_ca => 137.68179458376923,
sys.ltr_ca => 23.566070750638612,
sys.x_n1 => 0.02051602926129303,
sys.x_p3 => 0.010595307782827907,
sys.x_p2 => 0.012142111983813845,
sys.x_p1 => 0.00647030286225972,
sys.x_p0 => 0.007276352366632273,
sys.crp_ic => 16676.918901747533,
sys.crp_i => 16746.777783504043,
sys.adp_ic => 37.75097704519522,
sys.x_yca => 0.7245465013652926,
sys.cca4_lcc => 9.809958085343013e-23,
sys.cca3_lcc => 3.120039381551406e-17,
sys.cca2_lcc => 3.721593607615466e-12,
sys.cca1_lcc => 1.973012869511327e-7,
sys.cca0_lcc => 0.003922560427831359,
sys.o_lcc => 1.106771606240648e-23,
sys.c4_lcc => 9.725127301109256e-23,
sys.c3_lcc => 1.2377646305513033e-16,
sys.c2_lcc => 5.906074175046752e-11,
sys.c1_lcc => 1.2524994845840689e-5,
sys.pc2_ryr => 0.1268315869669829,
sys.po2_ryr => 2.8085358346208635e-9,
sys.po1_ryr => 0.00017254605527242085,
sys.O2 => 5.920151510560909,
sys.sox_m => 0.11424499555123306,
sys.dpsi => 171.47400536794677,
sys.nadh_m => 440.9173150524607,
sys.adp_m => 6.225787882031608,
sys.adp_i => 37.12263459699566,
sys.ca_m => 0.8518010873552284,
sys.ca_ss => 0.19792262061120025,
sys.ca_jsr => 1238.3749208770616,
sys.ca_nsr => 1242.3665362879096,
sys.ca_i => 0.19518567290699418,
sys.k_i => 144972.08735519781,
sys.na_i => 10213.437817331524,
sys.vm => -85.3671876586454,

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.3
Commit 966d0af0fdf (2025-12-15 11:20 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 12 × Intel(R) Core(TM) i7-8700 CPU @ 3.20GHz
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, skylake)
  GC: Built with stock GC
Threads: 2 default, 1 interactive, 2 GC (on 12 virtual cores)
Environment:
  JULIA_CI = true
  JULIA_CONDAPKG_OFFLINE = true
  LD_LIBRARY_PATH = /home/github/actions-runner-1/_work/_tool/Python/3.14.2/x64/lib
  JULIA_NUM_THREADS = 2
  JULIA_PROJECT = /home/github/actions-runner-1/_work/ecme-rirr-dox/ecme-rirr-dox/Project.toml
  JULIA_DEPOT_PATH = /home/github/.julia:/home/github/actions-runner-1/_work/_tool/julia/1.12.3/x64/local/share/julia:/home/github/actions-runner-1/_work/_tool/julia/1.12.3/x64/share/julia
  JULIA_CONDAPKG_BACKEND = Null
  JULIA_LOAD_PATH = @:@v#.#:@stdlib
using Pkg
Pkg.status()
Project ECMEDox v0.4.0
Status `~/actions-runner-1/_work/ecme-rirr-dox/ecme-rirr-dox/Project.toml`
  [459566f4] DiffEqCallbacks v4.11.0
  [0b91fe84] DisplayAs v0.1.6
  [f6369f11] ForwardDiff v1.3.0
  [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.12.0
  [1dea7af3] OrdinaryDiffEq v6.105.0
  [91a5bcdd] Plots v1.41.3
  [33c8b6b6] ProgressLogging v0.1.6
  [9672c7b4] SteadyStateDiffEq v2.8.0
⌅ [0c5d862f] Symbolics v6.58.0
  [37e2e46d] LinearAlgebra v1.12.0
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated`

This notebook was generated using Literate.jl.

Back to top