Hypoxia reperfusion simulation

Hypoxia reperfusion simulation#

using ProgressLogging
using OrdinaryDiffEq
using ModelingToolkit
using DiffEqCallbacks
using ECMEDox
using ECMEDox: second, Hz, μM, nM, build_stim_callbacks
using Plots
using DisplayAs: PNG
Plots.default(lw=1.5)

tend = 100.0second
bcl = 1.0second
@named sys = build_model(; bcl, tend)
u0 = build_u0(sys)
sts = unknowns(sys)
alg = KenCarp47()
@unpack O2 = sys
prob = ODEProblem(sys, [u0;  O2 => 6nM], tend)

reoxygen! = (integrator) -> begin
    integrator.ps[sys.O2] = 6μM
    set_proposed_dt!(integrator, 1e-6)
end

cbs = PresetTimeCallback(50.0second, reoxygen!)

┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/sc64L/src/variable.jl:587
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/sc64L/src/variable.jl:587
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/sc64L/src/variable.jl:587
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/sc64L/src/variable.jl:587
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/sc64L/src/variable.jl:587
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/sc64L/src/variable.jl:587
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/sc64L/src/variable.jl:587
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/sc64L/src/variable.jl:587
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/sc64L/src/variable.jl:587
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/sc64L/src/variable.jl:587
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/sc64L/src/variable.jl:587
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/sc64L/src/variable.jl:587
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/sc64L/src/variable.jl:587
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/sc64L/src/variable.jl:587
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/sc64L/src/variable.jl:587
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/sc64L/src/variable.jl:587
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/sc64L/src/variable.jl:587
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/sc64L/src/variable.jl:587
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/sc64L/src/variable.jl:587
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/sc64L/src/variable.jl:587
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/sc64L/src/variable.jl:587

SciMLBase.DiscreteCallback{DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Main.var"##230".var"#1#2"}, Main.var"##230".var"#1#2", DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Main.var"##230".var"#1#2"}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing}(DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Main.var"##230".var"#1#2"}([50000.0], true, SciMLBase.INITIALIZE_DEFAULT, Main.var"##230".var"#1#2"()), Main.var"##230".var"#1#2"(), DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Main.var"##230".var"#1#2"}([50000.0], true, SciMLBase.INITIALIZE_DEFAULT, Main.var"##230".var"#1#2"()), SciMLBase.FINALIZE_DEFAULT, Bool[1, 1], nothing)

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

  6.389516 seconds (21.16 M allocations: 1.020 GiB, 3.99% gc time, 82.68% compilation time)

retcode: Success
Interpolation: 3rd order Hermite
t: 10983-element Vector{Float64}:
      0.0
      1.0e-6
      1.1e-5
      0.00011099999999999999
      0.0011109999999999998
      0.0098657916786707
      0.02266582512295119
      0.04037540120694326
      0.06531031960839598
      0.10106523201078316
      ⋮
  99572.6383877471
  99620.61361158677
  99671.83555021304
  99728.58061877913
  99789.68385524656
  99853.4504871416
  99923.6668897242
 100000.0
 100000.0
u: 10983-element Vector{Vector{Float64}}:
 [0.3553843722148506, 162.59478245522047, 507.9034833706153, 14.892129267047379, 48.65695610700525, 0.9114480160471111, 0.19409331392839435, 1205.8091362833657, 1209.9687795381885, 144753.168650273  …  165.09011792191754, 91.87929617274456, 91.91725390039272, 1842.3823360869799, 131.3763898444155, 0.024427237952508167, 16.995660402659126, 1.9070224132868452, 0.0006449826277658302, 0.000865771566630073]
 [0.3553843443914563, 162.59480862690788, 507.90348338598517, 14.892129266274098, 48.65695610503596, 0.9114480160053736, 0.19409331391211282, 1205.8091361114257, 1209.9687793624942, 144753.16865524123  …  165.09013968705722, 91.8792795886062, 91.91725389999543, 1842.3823360868837, 131.37638984454472, 0.02446626519672757, 16.995660403658945, 1.907022413290002, 0.000644982627768534, 0.0008657715666339692]
 [0.35538406615764745, 162.59507019682923, 507.90348353969296, 14.892129258410554, 48.65695608533515, 0.9114480155879974, 0.19409331374930405, 1205.8091343920257, 1209.9687776055516, 144753.1687049235  …  165.09035727099683, 91.87911384816927, 91.91725387595574, 1842.38233605882, 131.37638984579925, 0.02485638711038522, 16.995660413817895, 1.9070224133215725, 0.0006449826277955726, 0.0008657715666850282]
 [0.3553812838330149, 162.59767125253344, 507.90348507773643, 14.892129166908138, 48.65695588753648, 0.9114480114141137, 0.1940933121308462, 1205.8091171980263, 1209.9687600361272, 144753.1692017371  …  165.09252638216864, 91.87746650761272, 91.91725163882411, 1842.3823330782172, 131.37638985462908, 0.028742592481314784, 16.99566053101768, 1.9070224136372766, 0.0006449826280660412, 0.0008657715684003225]
 [0.35535346192534245, 162.62226787374556, 507.9035006286395, 14.892127011583561, 48.65695383486602, 0.9114479696635892, 0.19409330030854405, 1205.8089452580962, 1209.9685843419768, 144753.17416897893  …  165.11356181599476, 91.8619692668177, 91.91703926770523, 1842.3820430034932, 131.37638958397153, 0.0661416025372887, 16.995662875593766, 1.9070224167943177, 0.000644982630851272, 0.0008657717011837188]
 [0.3551099865906897, 162.7557188386647, 507.90365997248773, 14.892037708390838, 48.656932478239895, 0.9114476034857547, 0.19409361484355522, 1205.8074399646596, 1209.967046200629, 144753.2175885664  …  165.2568262248503, 91.78471809215935, 91.90576248230506, 1842.364838820911, 131.37636634808047, 0.30235374646818314, 16.99569152370048, 1.9070224444337434, 0.0006449826994812201, 0.0008657791670530364]
 [0.35475430928565976, 162.813816826211, 507.9039375761985, 14.891798566906674, 48.65690119127446, 0.9114470671023076, 0.1940953038634002, 1205.8052391504918, 1209.9647975385544, 144753.28086542952  …  165.3853942978475, 91.77504854720662, 91.8826750675944, 1842.3185940917929, 131.37629731290062, 0.4675778437527647, 16.9957058669279, 1.9070224848467254, 0.0006449830884725159, 0.0008657983074115731]
 [0.3542627468960701, 162.81317112996342, 507.90433956060053, 14.891428867053323, 48.6568730450779, 0.9114463246437746, 0.19409985154400713, 1205.802194229248, 1209.961687124275, 144753.36805375997  …  165.4966893328008, 91.82340378464711, 91.86787393366978, 1842.2522705001788, 131.37618309543382, 0.5486304638324028, 16.995706509281423, 1.9070225407736152, 0.0006449842737589018, 0.000865823383255989]
 [0.3535716392289966, 162.78400468729492, 507.9048931473707, 14.890949309920874, 48.65687185942695, 0.9114452797208534, 0.19411043980802584, 1205.7979070689607, 1209.957310086573, 144753.4902216429  …  165.61357702171358, 91.90074777515314, 91.88343305984098, 1842.1752497175924, 131.37602901556588, 0.5760745213763646, 16.995700843305062, 1.907022619564239, 0.0006449868332325069, 0.0008658441990307148]
 [0.35258266227192986, 162.73850289401952, 507.90564518764376, 14.89038739378337, 48.65694663914088, 0.9114437827839661, 0.1941339972302191, 1205.7917597082092, 1209.9510412703776, 144753.66449940953  …  165.76412678555891, 91.9912611484821, 91.94240215786876, 1842.0871345283722, 131.37585670693738, 0.5816571652816934, 16.995692015984513, 1.9070227326643439, 0.0006449908306082892, 0.0008658403307081667]
 ⋮
 [0.05018021656021882, 141.61009958193648, 3.812733380881629, 196.06486597149075, 364.63303388221277, 1.031043743663059, 0.20106626799676774, 1162.7054984545184, 1172.0978703679123, 144377.70860671005  …  249.15703976632918, 12.127905888422976, 12.169439307307261, 1949.0644000451207, 77.348360982622, 0.1839210613326413, 16.96601638292397, 0.10091261032832981, 3.541355657606755e-5, 4.5690537971619385e-5]
 [0.05031683702537404, 141.6463975849261, 3.838872044184408, 195.14353376728928, 362.995286137327, 1.0294961110076737, 0.1961458431200156, 1157.1730509528948, 1165.452976937137, 144379.75320037565  …  249.02867756088767, 12.178291388190639, 12.219851063272083, 1948.9734624788746, 77.42893174957943, 0.1844045498037318, 16.966125050268356, 0.10093597046512974, 3.5527611869204596e-5, 4.583792565273334e-5]
 [0.050452166715873124, 141.68236978418344, 3.864869082783157, 194.2384959400264, 361.23079676654226, 1.0277750324703445, 0.19203138096671538, 1149.703089806434, 1157.033383333367, 144382.03798860597  …  248.90173461510787, 12.228140325181679, 12.269725767769042, 1948.883429044269, 77.50877290947943, 0.18488363238511707, 16.966232105828038, 0.10096625689872968, 3.564081939836078e-5, 4.598422466546161e-5]
 [0.050590549014511946, 141.71946392111448, 3.891775032116921, 193.31133827047174, 359.31858038340715, 1.0258043074305951, 0.18855951884909117, 1140.2120105966671, 1146.72300536927, 144384.66503288294  …  248.77238412398808, 12.27904305567628, 12.32065410608626, 1948.7915959287086, 77.59009534640772, 0.18537362257083106, 16.9663418292416, 0.10100564068818238, 3.5757112494158015e-5, 4.613451770802935e-5]
 [0.050727414548077554, 141.7566512347493, 3.9188548139110915, 192.38670335934418, 357.3471394905464, 1.0236259831080772, 0.18576947094521004, 1129.2129048006577, 1135.0509789546045, 144387.5798888271  …  248.64506586786615, 12.329311196861548, 12.370946530720103, 1948.7011085317954, 77.67000339714447, 0.18585829936763004, 16.96645115998419, 0.10105420995575357, 3.587285811876924e-5, 4.628410923398774e-5]
 [0.05085833844488634, 141.7927899556108, 3.9452792339527916, 191.491953626315, 355.4010718339764, 1.0213073789588678, 0.18359288153132264, 1117.4182664461518, 1122.7216559208882, 144390.69389203624  …  248.52393312459628, 12.377322855452462, 12.41898028376393, 1948.6149201569835, 77.74585154046055, 0.1863219627590486, 16.96655678768653, 0.10111097748352228, 3.598436489261187e-5, 4.642822769488308e-5]
 [0.05099011829347106, 141.82977833203236, 3.9724379550365145, 190.57945399175054, 353.3865497320151, 1.0187156548432792, 0.18177508619335633, 1104.506504526632, 1109.3579001527924, 144394.18801504045  …  248.40270424340798, 12.425573609071753, 12.467252059636031, 1948.5285604446121, 77.8215631627396, 0.1867886737802605, 16.966664277295756, 0.10117984963225125, 3.6097439543581664e-5, 4.6574377772930564e-5]
 [0.051120496132547, 141.86702461009494, 3.999897843590313, 189.66346410197224, 351.3325748727978, 1.015864968953399, 0.1802393566203331, 1090.9000650824053, 1095.3684139990175, 144398.04728139774  …  248.28348573956222, 12.47323314328211, 12.514931104174076, 1948.4435234277237, 77.89581965010498, 0.1872504382166418, 16.966771873456967, 0.10126135672512203, 3.621019206018013e-5, 4.672011691289023e-5]
 [0.051120496132547, 141.86702461009494, 3.999897843590313, 189.66346410197224, 351.3325748727978, 1.015864968953399, 0.1802393566203331, 1090.9000650824053, 1095.3684139990175, 144398.04728139774  …  248.28348573956222, 12.47323314328211, 12.514931104174076, 1948.4435234277237, 77.89581965010498, 0.1872504382166418, 16.966771873456967, 0.10126135672512203, 3.621019206018013e-5, 4.672011691289023e-5]

@unpack vm, dpsi, atp_i, adp_i = sys
pl_mmp = plot(sol, idxs=dpsi, lab=false, title="(A) Mito. memb. potential", xlabel="", ylabel="Voltage (mV)")
pl_mmp = vline!(pl_mmp, [50.0second], lines=(:dash, :black), lab=false)
pl_vm = plot(sol, idxs=vm, lab=false, title="(B) Action potential", xlabel="", ylabel="Voltage (mV)")
pl_vm = vline!(pl_vm, [50.0second], lines=(:dash, :black), lab=false)
pl_atp = plot(sol, idxs=atp_i/1000, lab=false, title="(C) ATP", xlabel="", ylabel="Conc. (mM)")
pl_atp = vline!(pl_atp, [50.0second], lines=(:dash, :black), lab=false)
@unpack cit, isoc, oaa, akg, scoa, suc, fum, mal = sys
pl_cac = plot(sol, idxs=[cit, isoc, oaa, akg, scoa, suc, fum, mal], legend=:right, title="(D) CAC intermediates", xlabel="Time (ms)", ylabel="Conc. (μM)")
pl_cac = vline!(pl_cac, [50.0second], lines=(:dash, :black), lab=false)
@unpack Q_n, SQn, QH2_n, QH2_p, Q_p, fes_ox, fes_rd, cytc_ox, cytc_rd = sys
pl_q = plot(sol, idxs=[Q_n, Q_p, SQn, QH2_n, QH2_p], title="(E) Q cycle", legend=:left, xlabel="Time (ms)", ylabel="Conc. (μM)")
pl_q = vline!(pl_q, [50.0second], lines=(:dash, :black), lab=false)
pl_ros = plot(sol, idxs=100 * sys.vROS / (sys.vO2 + sys.vROS), title="(F) ROS generation", lab=false, xlabel="Time (ms)", ylabel="Fraction of O2 consumption (%)")
pl_ros = vline!(pl_ros, [50.0second], lines=(:dash, :black), lab=false)
plot(pl_mmp, pl_vm, pl_atp, pl_cac, pl_q, pl_ros, size=(1200, 800)) |> PNG

_images/af4f3df40dfc524437bac563469cf260e411647895b4f85287325ecd0fb65776.png

plot(sol, idxs=sys.nadh_m, title="NADH (mito)")

_images/6c55fcd1c50bd236b655cb534cab48b9803f1aa3387aecdb3c3f098ad94fb212.png

plot(sol, idxs=[sys.vROSC1, sys.vROSIf, sys.vROSIq, sys.vROSC3], ylims=(0.0, 0.05))

_images/bfc0fbfd73dc4eeb7b38cdb1965eb8ae76354480ed488104576efd50edaa5898.png

plot(sol, idxs=[sys.sox_m, sys.sox_i])

_images/0547a10f6e1b8e6505a41806bc9e59d9008eecfa7fc7028245cdeba55aae5a8d.png

plot(sol, idxs=[sys.vSDH])

_images/5dad0af8ad36f1fd90f7deb0e1a24651f4e6f82f57f599d7803fbd19f793faf7.png

This notebook was generated using Literate.jl.