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"#2#4"{Float64, Float64, Symbolics.Num}}, ECMEDox.var"#2#4"{Float64, Float64, Symbolics.Num}, DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), ECMEDox.var"#2#4"{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"#3#5"{Float64, Float64, Symbolics.Num}}, ECMEDox.var"#3#5"{Float64, Float64, Symbolics.Num}, DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), ECMEDox.var"#3#5"{Float64, Float64, Symbolics.Num}}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}, 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, Tuple{}}}}((), (SciMLBase.DiscreteCallback{DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), ECMEDox.var"#2#4"{Float64, Float64, Symbolics.Num}}, ECMEDox.var"#2#4"{Float64, Float64, Symbolics.Num}, DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), ECMEDox.var"#2#4"{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"#2#4"{Float64, Float64, Symbolics.Num}}(0.0:1000.0:100000.0, true, SciMLBase.INITIALIZE_DEFAULT, ECMEDox.var"#2#4"{Float64, Float64, Symbolics.Num}(-80.0, 0.5, iStim(t))), ECMEDox.var"#2#4"{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"#2#4"{Float64, Float64, Symbolics.Num}}(0.0:1000.0:100000.0, true, SciMLBase.INITIALIZE_DEFAULT, ECMEDox.var"#2#4"{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"#3#5"{Float64, Float64, Symbolics.Num}}, ECMEDox.var"#3#5"{Float64, Float64, Symbolics.Num}, DiffEqCallbacks.PresetTimeFunction{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, typeof(SciMLBase.INITIALIZE_DEFAULT), ECMEDox.var"#3#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"#3#5"{Float64, Float64, Symbolics.Num}}(0.5:1000.0:100000.5, true, SciMLBase.INITIALIZE_DEFAULT, ECMEDox.var"#3#5"{Float64, Float64, Symbolics.Num}(0.0, 0.5, iStim(t))), ECMEDox.var"#3#5"{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"#3#5"{Float64, Float64, Symbolics.Num}}(0.5:1000.0:100000.5, true, SciMLBase.INITIALIZE_DEFAULT, ECMEDox.var"#3#5"{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"##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, Tuple{}}(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)
10.241362 seconds (17.58 M allocations: 858.493 MiB, 2.09% gc time, 81.33% 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.22496703875
99637.07314145607
99689.06897762191
99746.66745754494
99808.71359964785
99873.7591540077
99945.74659208026
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.20648429974023197, 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.2064647957809869, 1295.525001638845, 1299.782408146414, 0.20379328553620607, 145551.4924783092, 10276.491167946646, -79.19728404284417]
⋮
[0.12027092431973568, 4.223823539738923e-5, 5.451968684402189e-5, 245.79952985957755, 281.40906387116934, 197.34055251186768, 24.567609750742918, 42.45540792098763, 257.60432514800436, 0.1349757637258448 … 73.45008976414687, 119.95957838955684, 0.9478534741182, 0.1980199436067818, 1213.2803939941657, 1220.1555581733785, 0.19416461894577172, 144971.53157212166, 9819.848521471105, -85.22261833308451]
[0.12030206562381095, 4.240163163980496e-5, 5.4731127496092556e-5, 245.6365173689221, 281.3067272656202, 197.13491046104446, 24.583899641463663, 42.531760195821356, 257.51064177159657, 0.1356073944014109 … 72.9443530787421, 119.58829731011245, 0.9463112605427088, 0.19347032513523468, 1204.440247123602, 1210.317297244523, 0.18978252902990944, 144973.46417192923, 9827.837635295988, -85.329879731914]
[0.12034325456671006, 4.257229536195235e-5, 5.4951966752481524e-5, 245.47880807089714, 281.20766716597893, 196.93612557409492, 24.599604671638886, 42.60583539982524, 257.4198504059755, 0.13622193316325115 … 72.45835416553639, 119.20101082216335, 0.9446018294171918, 0.18970160137687025, 1193.6960324027991, 1198.7442468534832, 0.186182177385751, 144975.61704231307, 9836.004504362094, -85.41527973733105]
[0.12039815286935102, 4.2756404702653756e-5, 5.519019957020509e-5, 245.32037585732502, 281.10811457025045, 196.7367646871714, 24.615344383577515, 42.68090512864395, 257.3280177215773, 0.13684450955806987 … 71.96846618965749, 118.78921068798502, 0.9426454288433482, 0.18652458226278512, 1180.8507792944533, 1185.2029745508798, 0.18318054059636002, 144978.09173845281, 9844.738182434578, -85.48531234610535]
[0.12046762128499616, 4.294798286547153e-5, 5.54380935522815e-5, 245.1656138039357, 281.01083697615985, 196.54242052338304, 24.630729174277086, 42.75513080933662, 257.23739758352923, 0.1374587413590896 … 71.4863792874316, 118.36921676750627, 0.9404827492310588, 0.1839544174163672, 1166.5273275775992, 1170.3245308763524, 0.18078718253625772, 144980.83861634156, 9853.8671656694, -85.54096317430904]
[0.12055122953255418, 4.314063909605664e-5, 5.56873816562343e-5, 245.01791282114925, 280.91797039636083, 196.35733293520096, 24.645446690957055, 42.82687348763645, 257.1499633492702, 0.13805081944050254 … 71.02273721591624, 117.95483513534204, 0.9381702425492109, 0.18191236517604803, 1151.4792136035553, 1154.8478944883414, 0.17891804697027924, 144983.78760429577, 9863.202697737273, -85.5847706325535]
[0.12065567984554564, 4.334374523810341e-5, 5.5950193157766994e-5, 244.86868363674026, 280.82411584599265, 196.17072804488808, 24.660331776099675, 42.90019636966983, 257.06076174480205, 0.13865452449025542 … 70.55101111201543, 117.52392684748811, 0.9355711847057938, 0.1801632017161251, 1135.1946502746507, 1138.2125806145696, 0.17734664627635835, 144987.11589524857, 9873.32001182804, -85.62222076018885]
[0.12074205045371005, 4.348981912984151e-5, 5.613920963544078e-5, 244.7645240051903, 280.7585909740003, 196.04073115685142, 24.670691766667577, 42.95179374892503, 256.9981065448617, 0.1390787847481167 … 70.21981233180205, 117.21564442565942, 0.9335912342088238, 0.17907988168334235, 1123.3374022340115, 1126.1492912806702, 0.17638859966887016, 144989.66106902852, 9880.824977115573, -85.6454914224471]
[0.12074205045371005, 4.348981912984151e-5, 5.613920963544078e-5, 244.7645240051903, 280.7585909740003, 196.04073115685142, 24.670691766667577, 42.95179374892503, 256.9981065448617, 0.1390787847481167 … 70.21981233180205, 117.21564442565942, 0.9335912342088238, 0.17907988168334235, 1123.3374022340115, 1126.1492912806702, 0.17638859966887016, 144989.66106902852, 9880.824977115573, -85.6454914224471]
@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.