Hypoxia reperfusion simulation#
using ProgressLogging
using OrdinaryDiffEq
using ModelingToolkit
using DiffEqCallbacks
using ECMEDox
using ECMEDox: second, Hz, μM, nM, mV
using Plots
using DisplayAs: PNG
Plots.default(lw=1.5)
tend = 100.0second
bcl = 1.0second
@named sys = build_model()
u0 = build_u0(sys)
sts = unknowns(sys)
alg = KenCarp47()
@unpack iStim = sys
stim! = build_stim_callbacks(iStim, tend; period=bcl)
prob = ODEProblem(sys, [u0; sys.O2_o => 6nM], tend)
reoxygen! = (integrator) -> begin
integrator.ps[sys.O2_o] = 6μM
set_proposed_dt!(integrator, 1e-6)
end
cbs = CallbackSet(stim!, PresetTimeCallback(50.0second, reoxygen!))
SciMLBase.CallbackSet{Tuple{}, Tuple{SciMLBase.DiscreteCallback{DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), ECMEDox.var"#4#5"{Float64, Float64, Symbolics.Num}}, ECMEDox.var"#4#5"{Float64, Float64, Symbolics.Num}, DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), ECMEDox.var"#4#5"{Float64, Float64, Symbolics.Num}}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}, SciMLBase.DiscreteCallback{DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), ECMEDox.var"#6#7"{Float64, Float64, Symbolics.Num}}, ECMEDox.var"#6#7"{Float64, Float64, Symbolics.Num}, DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), ECMEDox.var"#6#7"{Float64, Float64, Symbolics.Num}}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}, SciMLBase.DiscreteCallback{DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Main.var"##295".var"#2#3"}, Main.var"##295".var"#2#3", DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Main.var"##295".var"#2#3"}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}}}((), (SciMLBase.DiscreteCallback{DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), ECMEDox.var"#4#5"{Float64, Float64, Symbolics.Num}}, ECMEDox.var"#4#5"{Float64, Float64, Symbolics.Num}, DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), ECMEDox.var"#4#5"{Float64, Float64, Symbolics.Num}}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}(DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), ECMEDox.var"#4#5"{Float64, Float64, Symbolics.Num}}(0.0:1000.0:100000.0, true, SciMLBase.INITIALIZE_DEFAULT, ECMEDox.var"#4#5"{Float64, Float64, Symbolics.Num}(-80.0, 0.5, iStim(t))), ECMEDox.var"#4#5"{Float64, Float64, Symbolics.Num}(-80.0, 0.5, iStim(t)), DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), ECMEDox.var"#4#5"{Float64, Float64, Symbolics.Num}}(0.0:1000.0:100000.0, true, SciMLBase.INITIALIZE_DEFAULT, ECMEDox.var"#4#5"{Float64, Float64, Symbolics.Num}(-80.0, 0.5, iStim(t))), SciMLBase.FINALIZE_DEFAULT, Bool[1, 1], nothing, ()), SciMLBase.DiscreteCallback{DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), ECMEDox.var"#6#7"{Float64, Float64, Symbolics.Num}}, ECMEDox.var"#6#7"{Float64, Float64, Symbolics.Num}, DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), ECMEDox.var"#6#7"{Float64, Float64, Symbolics.Num}}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}(DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), ECMEDox.var"#6#7"{Float64, Float64, Symbolics.Num}}(0.5:1000.0:100000.5, true, SciMLBase.INITIALIZE_DEFAULT, ECMEDox.var"#6#7"{Float64, Float64, Symbolics.Num}(0.0, 0.5, iStim(t))), ECMEDox.var"#6#7"{Float64, Float64, Symbolics.Num}(0.0, 0.5, iStim(t)), DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), ECMEDox.var"#6#7"{Float64, Float64, Symbolics.Num}}(0.5:1000.0:100000.5, true, SciMLBase.INITIALIZE_DEFAULT, ECMEDox.var"#6#7"{Float64, Float64, Symbolics.Num}(0.0, 0.5, iStim(t))), SciMLBase.FINALIZE_DEFAULT, Bool[1, 1], nothing, ()), SciMLBase.DiscreteCallback{DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Main.var"##295".var"#2#3"}, Main.var"##295".var"#2#3", DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Main.var"##295".var"#2#3"}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}(DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Main.var"##295".var"#2#3"}([50000.0], true, SciMLBase.INITIALIZE_DEFAULT, Main.var"##295".var"#2#3"()), Main.var"##295".var"#2#3"(), DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), Main.var"##295".var"#2#3"}([50000.0], true, SciMLBase.INITIALIZE_DEFAULT, Main.var"##295".var"#2#3"()), 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)
12.621327 seconds (21.68 M allocations: 1.280 GiB, 1.99% gc time, 81.81% compilation time)
retcode: Success
Interpolation: 3rd order Hermite
t: 9127-element Vector{Float64}:
0.0
1.0e-6
1.1e-5
0.00011099999999999999
0.0011109999999999998
0.0082007461100045
0.018186650646576062
0.031947957368564726
0.05075735944316681
0.07717716533731177
⋮
99588.22521616989
99637.07341600092
99689.06927357841
99746.66778167286
99808.71394833988
99873.75953345661
99945.74700556856
100000.0
100000.0
u: 9127-element Vector{Vector{Float64}}:
[2.196435053752864, 0.0007273589705197244, 0.000981691808598968, 245.22202099807583, 281.0719264404058, 197.06621753790645, 34.41477140951692, 60.37913837798986, 229.2919970960533, 22.321405235243013 … 1.7002259496890006, 23.148552252520883, 0.6928662120941113, 0.20648605639609985, 1295.5364065086924, 1299.7940884726063, 0.20381438659933854, 145551.11317854663, 10276.480036234318, -85.29107641956527]
[2.1964350537444246, 0.0007273589704930469, 0.0009816918085575927, 245.2220201608852, 281.0719264411974, 197.06571474064046, 34.413522081818975, 60.380875446623286, 229.2915014805287, 22.320902445095992 … 1.7002259496141854, 23.148552251328372, 0.6928662120479763, 0.2064860563812371, 1295.536406360914, 1299.7940883214064, 0.20381438658633685, 145551.11318351026, 10276.48003640626, -85.29099641994007]
[2.196435053660033, 0.0007273589702262693, 0.0009816918077433308, 245.22201178911322, 281.0719264463143, 197.06068746841765, 34.401035657524744, 60.398238636137705, 229.28654599490577, 22.315875848604467 … 1.7002259489873959, 23.148552239710824, 0.6928662115866655, 0.2064860562325451, 1295.5364048831284, 1299.794086809408, 0.2038143864524468, 145551.11323314617, 10276.480038125655, -85.29019642545443]
[2.196435052816117, 0.0007273589675557858, 0.000981691759697, 245.22192808483354, 281.0719262178529, 197.01048468391616, 34.27685412851235, 60.57112348934708, 229.23705800590776, 22.265740120278917 … 1.700225954954034, 23.14855215416808, 0.6928662069775376, 0.20648605467924058, 1295.5363901052715, 1299.794071689427, 0.2038143847280742, 145551.11372949643, 10276.480055315615, -85.28219665718146]
[2.1964350443769423, 0.0007273589382046185, 0.0009816874319058665, 245.22109235353753, 281.0718962429053, 196.51534552184557, 33.10083377880719, 62.22782944610671, 228.7487703255649, 21.777148664722105 … 1.7002273272600748, 23.1485542416841, 0.6928661612678696, 0.20648601087336682, 1295.5362423267495, 1299.7939204894888, 0.203814330723474, 145551.11869213346, 10276.48022681667, -85.202216589277]
[2.1964349845427154, 0.0007273578428589997, 0.0009815015326164894, 245.2152231707679, 281.0703366884667, 193.3265342469898, 27.42237585380103, 71.01466301494123, 225.5963664338579, 18.889006971615018 … 1.7003209983959373, 23.148685240598265, 0.6928658521264267, 0.20648429974023194, 1295.535194617031, 1299.792848464894, 0.203812573311686, 145551.1538318486, 10276.481421950844, -84.6360800777001]
[2.1964349002294767, 0.0007273504709519492, 0.000980977105944279, 245.20704847031118, 281.0646671150038, 189.58665395991727, 24.03777466558498, 78.06554454667612, 221.88546777241575, 16.06487462896432 … 1.700673404295367, 23.149058461025486, 0.6928654408072512, 0.2064799534029841, 1295.5337189257127, 1299.7913381709932, 0.20380822314203575, 145551.20320141272, 10276.483043682256, -83.84125168021893]
[2.196434783823074, 0.0007273244385408795, 0.0009800788277491348, 245.1957778624932, 281.05149164059424, 185.4084913756189, 22.891954527192137, 83.32395546756774, 217.7264422147432, 13.550596291918247 … 1.7013525027047107, 23.149687632585046, 0.6928648892732475, 0.20647369405881294, 1295.5316853294507, 1299.7892560337336, 0.20380200154835443, 145551.27101776408, 10276.485160911678, -82.75051592544995]
[2.196434623862003, 0.000727256848993229, 0.000978809822655273, 245.18008319472884, 281.02564745662505, 180.86166217210754, 23.169137298437885, 87.53576518973242, 213.1865635638576, 11.373425271240329 … 1.702318652533269, 23.150554649823363, 0.6928641384996909, 0.20646746294247903, 1295.5289057709726, 1299.786408757247, 0.2037958438859671, 145551.3633526615, 10276.487835521653, -81.26743606933216]
[2.1964343963246744, 0.0007271034468491387, 0.0009771117655533953, 245.15704996925584, 280.97789543613624, 175.90331302507056, 24.04333191114688, 91.57319134413137, 208.21641942771203, 9.347618887913466 … 1.7035430379792826, 23.151672306676424, 0.6928630755349369, 0.20646479578098687, 1295.525001638845, 1299.782408146414, 0.20379328553620607, 145551.4924783092, 10276.491167946646, -79.19728404284417]
⋮
[0.12027092445183382, 4.2238236232749e-5, 5.451968792502189e-5, 245.79952899119098, 281.40906332618795, 197.34055141630134, 24.567609837728103, 42.455408328127895, 257.6043246483349, 0.13497576708467654 … 73.45008704931259, 119.9595765077744, 0.9478534664713609, 0.19801991754163167, 1213.280353067846, 1220.1555115418853, 0.19416459375815703, 144971.5315817834, 9819.848562929035, -85.22261896243819]
[0.12030206581199274, 4.240163254667686e-5, 5.473112866960545e-5, 245.63651649795176, 281.3067267186652, 197.13490936241996, 24.583899728433085, 42.531760603795725, 257.5106412710814, 0.1356073977840336 … 72.94435039512851, 119.58829524963198, 0.9463112517209801, 0.19347030253910827, 1204.4401935815022, 1210.3172387353036, 0.1897825073482159, 144973.46418311982, 9827.837679059781, -85.32988025257954]
[0.12034325481639413, 4.257229631809589e-5, 5.4951967989715866e-5, 245.47880721808383, 281.20766663019094, 196.93612449985875, 24.599604756540344, 42.60583580174746, 257.419849913666, 0.1362219364986818 … 72.45835153896915, 119.20100866914255, 0.9446018195537365, 0.1897015826153815, 1193.6959685162778, 1198.7441788432502, 0.18618215955377138, 144975.61705486948, 9836.004549809282, -85.41528015571718]
[0.12039815319698156, 4.2756405718499516e-5, 5.5190200884673064e-5, 245.32037501127357, 281.10811403852733, 196.7367636234826, 24.61534446769529, 42.680905531806985, 257.32801722880527, 0.13684451289909288 … 71.96846356739741, 118.78920844247833, 0.9426454177151734, 0.18652456679716864, 1180.8507053914088, 1185.2028972778053, 0.18318052608542587, 144978.0917526649, 9844.738230610074, -85.48531268331296]
[0.12046762169639968, 4.2947983917213255e-5, 5.5438094913184036e-5, 245.1656129774874, 281.0108364565992, 196.5424194865609, 24.63072925665065, 42.75513120828692, 257.23739709678176, 0.1374587446559832 … 71.48637670554344, 118.36921449002136, 0.940482736980892, 0.18395440491831871, 1166.5272466986928, 1170.3244473318077, 0.18078717099748662, 144980.8386320388, 9853.867216100898, -85.54096344322149]
[0.12055123004300254, 4.314064019114079e-5, 5.56873830732198e-5, 245.01791200008248, 280.9179698800352, 196.35733190722573, 24.64544677301122, 42.82687388899221, 257.14996286041395, 0.1380508227485408 … 71.02273463210466, 117.95483280411915, 0.938170228983931, 0.1819123547442934, 1151.4791265361962, 1154.8478052759065, 0.17891803751131968, 144983.78762173615, 9863.202751413111, -85.58477085582035]
[0.12065568047148628, 4.3343746370588945e-5, 5.595019462316841e-5, 244.86868281961654, 280.8241153319974, 196.1707270241205, 24.660331857727176, 42.90019677333917, 257.0607612540378, 0.13865452781038276 … 70.55100852326352, 117.52392446059521, 0.9355711697137427, 0.18016319281062346, 1135.194558473665, 1138.2124870803175, 0.17734663835247194, 144987.11591459572, 9873.3200691929, -85.62222095105398]
[0.12074205044561526, 4.348981912603918e-5, 5.6139209630518225e-5, 244.7645240076186, 280.7585909755218, 196.0407311598304, 24.670691766500674, 42.95179374758633, 256.9981065463834, 0.13907878473741914 … 70.21981234222699, 117.21564443878461, 0.9335912342445001, 0.17907988168790495, 1123.3374023024055, 1126.1492913499235, 0.17638859967268336, 144989.66106908952, 9880.824976915428, -85.64549142230878]
[0.12074205044561526, 4.348981912603918e-5, 5.6139209630518225e-5, 244.7645240076186, 280.7585909755218, 196.0407311598304, 24.670691766500674, 42.95179374758633, 256.9981065463834, 0.13907878473741914 … 70.21981234222699, 117.21564443878461, 0.9335912342445001, 0.17907988168790495, 1123.3374023024055, 1126.1492913499235, 0.17638859967268336, 144989.66106908952, 9880.824976915428, -85.64549142230878]
@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 (%)", ylims=(0, 100))
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
MMP
@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)")
vline!(pl_mmp, [50.0second], lines=(:dash, :black), lab=false)
pl_mmp |> PNG
Action potential
pl_vm = plot(sol, idxs=vm, lab=false, title="(B) Action potential", xlabel="", ylabel="Voltage (mV)")
vline!(pl_vm, [50.0second], lines=(:dash, :black), lab=false)
pl_vm |> PNG
ATP
pl_atp = plot(sol, idxs=atp_i/1000, lab=false, title="(C) ATP", xlabel="", ylabel="Conc. (mM)")
vline!(pl_atp, [50.0second], lines=(:dash, :black), lab=false)
pl_atp |> PNG
CAC
@unpack cit, isoc, oaa, akg, scoa, suc, fum, mal, nadh_m = sys
pl_cac = plot(sol, idxs=[cit, isoc, oaa, akg, scoa, suc, fum, mal, nadh_m], legend=:right, title="(D) CAC intermediates", xlabel="Time (ms)", ylabel="Conc. (μM)")
vline!(pl_cac, [50.0second], lines=(:dash, :black), lab=false)
pl_cac |> PNG
CAC flux
pl = plot(sol, idxs=[sys.vSDH, sys.vAAT, sys.vMDH, sys.vIDH])
vline!(pl, [50.0second], lines=(:dash, :black), lab=false)
pl |> PNG
Q cycle
@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)")
vline!(pl_q, [50.0second], lines=(:dash, :black), lab=false)
pl_q |> PNG
Downstream of Q cycle
pl = plot(sol, idxs=[fes_ox, fes_rd, cytc_ox, cytc_rd], title="Q cycle downstream", legend=:left, xlabel="Time (ms)", ylabel="Conc. (μM)")
vline!(pl, [50.0second], lines=(:dash, :black), lab=false)
pl |> PNG
Proton pumping
pl = plot(sol, idxs=[sys.vHresC1, sys.vHresC3, sys.vHresC4], title="Proton pumping", legend=:left, xlabel="Time (ms)", ylabel="Rate (μM/ms)")
vline!(pl, [50.0second], lines=(:dash, :black), lab=false)
pl |> PNG
Mito oxygen concentration
pl = plot(sol, idxs=sys.O2)
vline!(pl, [50.0second], lines=(:dash, :black), lab=false)
pl |> PNG
O2 shunt
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 (%)", ylims=(0, 100))
vline!(pl_ros, [50.0second], lines=(:dash, :black), lab=false)
pl_ros |> PNG
Superoxide production
pl = plot(sol, idxs=[sys.vROSC1, sys.vROSIf, sys.vROSIq, sys.vROSC3], ylims=(0.0, 0.05))
vline!(pl, [50.0second], lines=(:dash, :black), lab=false)
pl |> PNG
Superoxide concentrations
pl = plot(sol, idxs=[sys.sox_m, sys.sox_i])
vline!(pl, [50.0second], lines=(:dash, :black), lab=false)
pl |> PNG
Semiquinone from Complex III Qo
pl = plot(sol, idxs=[sys.SQp])
vline!(pl, [50.0second], lines=(:dash, :black), lab=false)
pl |> PNG
This notebook was generated using Literate.jl.