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}, 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}, 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}}}((), (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}(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}(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}(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)
7.896823 seconds (17.55 M allocations: 870.389 MiB, 2.98% gc time, 71.05% compilation time)
retcode: Success
Interpolation: 3rd order Hermite
t: 11965-element Vector{Float64}:
0.0
1.0e-6
1.1e-5
0.00011099999999999999
0.0011109999999999998
0.011110999999999996
0.04420707379724256
0.12080260736933432
0.19431330755515736
0.33607770270628173
⋮
99614.8657235206
99665.10174983881
99721.42415050876
99782.06786449286
99845.6883962849
99915.6369614186
99991.89505097178
100000.0
100000.0
u: 11965-element Vector{Vector{Float64}}:
[1.4355864646916074, 0.0004949419698468167, 0.0006578256307800626, 209.31680906407098, 257.044146799667, 156.82453780216784, 54.88374367717081, 128.68082835522182, 136.13456850129444, 168.9705244901982 … 11.000668810590449, 43.96538500887081, 0.8926956182412644, 0.19547857075761974, 1217.9822251524322, 1222.0613051409616, 0.19267665443392065, 144821.2542377956, 10200.108762034, -85.39304888131878]
[1.4355864646996925, 0.0004949419698519347, 0.0006578256307870891, 209.31680906451868, 257.0441468000075, 156.8245378028954, 54.88374367671794, 128.68082835499527, 136.13456850205296, 168.97052448772078 … 11.000668918503342, 43.96538500707363, 0.892695618199064, 0.1954785707415604, 1217.982224981597, 1222.0613049664905, 0.19267665442002282, 144821.2542427633, 10200.108762189264, -85.3929688816996]
[1.4355864647805472, 0.0004949419699031138, 0.0006578256308573539, 209.31680906898904, 257.0441468034126, 156.82453781017085, 54.88374367218921, 128.6808283527297, 136.13456850963792, 168.97052446294666 … 11.000669997631693, 43.96538498910568, 0.8926956177770586, 0.19547857058098675, 1217.9822232732456, 1222.061303221779, 0.19267665428149633, 144821.25429244002, 10200.108763741864, -85.39216888730721]
[1.4355864655890933, 0.0004949419704149057, 0.0006578256315599976, 209.31680911301157, 257.0441468374625, 156.82453788292526, 54.88374362690356, 128.6808283300725, 136.13456858548733, 168.97052421520564 … 11.000680788856199, 43.96538480981478, 0.8926956135570056, 0.19547856898400512, 1217.9822061897326, 1222.0612857746653, 0.1926766529414278, 144821.2547891986, 10200.108779264028, -85.3841691232638]
[1.4355864736745558, 0.0004949419755328263, 0.0006578256385829357, 209.3168094851509, 257.0441471779335, 156.82453861034912, 54.88374317566813, 128.68082810198914, 136.1345693438671, 168.9705217378717 … 11.000788695204921, 43.96538305572277, 0.8926955713564741, 0.19547855647577966, 1217.9820353546659, 1222.0611113036139, 0.1926766440593151, 144821.2597558932, 10200.108934100916, -85.30418942657788]
[1.4355865545293158, 0.000494942026695752, 0.000657825705333404, 209.31680640647437, 257.04415055503256, 156.82454578043644, 54.88373997078312, 128.68082460812602, 136.13457682902538, 168.97049701884168 … 11.001867171949922, 43.96536935537744, 0.8926951493517562, 0.19547888026211724, 1217.9803270101745, 1222.0593666221212, 0.19267700604969926, 144821.30933639948, 10200.110444052265, -84.50614320939705]
[1.435586822128531, 0.0004949421920598613, 0.0006578257233988584, 209.31670844322647, 257.0441601175139, 156.82456539433733, 54.883772725979895, 128.68077336447666, 136.1345977513635, 168.9704162650728 … 11.005429030706633, 43.96537173485913, 0.8926937527269715, 0.1954857462736804, 1217.9746731407804, 1222.053593815092, 0.1926839973076688, 144821.47243726387, 10200.114946272264, -81.88548232853988]
[1.4355874414489875, 0.000494942375547163, 0.0006578221497308882, 209.31596977101793, 257.04415364967997, 156.82457228152057, 54.884182403767745, 128.680354045614, 136.13461138457174, 168.97022571611163 … 11.013630272136762, 43.965626969826424, 0.892690521316143, 0.19553500284837125, 1217.9615886453198, 1222.0402579837494, 0.19273353116101058, 144821.84595357705, 10200.122495467709, -75.91031344424186]
[1.435588035796651, 0.0004949416133115159, 0.0006578093784375139, 209.31459525146474, 257.0440745442874, 156.82451893145992, 54.885033453158705, 128.67954224910852, 136.13457328749223, 168.97002676846574 … 11.021449711947763, 43.96614229808763, 0.8926874223162142, 0.1956247362731244, 1217.949031755149, 1222.0275228575258, 0.19282350366355466, 144822.20241835254, 10200.126062174759, -70.24094585579968]
[1.4355891814057113, 0.0004949327047120863, 0.000657740501143926, 209.31011177697363, 257.04358376648804, 156.82421800466952, 54.887984939112656, 128.67681214956957, 136.1343448737221, 168.969579805377 … 11.036403649622217, 43.967665601728896, 0.8926814570798979, 0.1959102078619084, 1217.9248178480761, 1222.003272840465, 0.1931093185020087, 144822.8932537627, 10200.123727809092, -59.327549973580574]
⋮
[0.11595018930219796, 4.081952703175054e-5, 5.2683392437673446e-5, 255.5221841189699, 287.3623947562443, 209.27737937635814, 31.296225545905926, 48.90180855131449, 244.17017195592976, 14.274311796991874 … 145.33338682521926, 236.44979323323747, 1.0234474792761186, 0.19810483011875207, 1180.0515867034458, 1188.1827058975239, 0.1939191217389556, 144448.86568853448, 10037.651344953763, -85.23908361424249]
[0.1159852101603092, 4.0951414230434175e-5, 5.285396955556529e-5, 255.43454424335127, 287.3092098074619, 209.1602021498085, 31.331510804657704, 48.99101060301585, 244.04350016555733, 14.330496114391906 … 144.63709923029126, 235.52783852951944, 1.0217531502888175, 0.194005021228751, 1172.3433009778048, 1179.5294979396117, 0.18999658563871225, 144451.03477835597, 10045.574544201323, -85.33196249792037]
[0.11603133440305753, 4.1091022563579975e-5, 5.303454183485701e-5, 255.3423140045327, 287.2532137118819, 209.03677821031818, 31.36908751087847, 49.08554885563941, 243.90906046490144, 14.390139149252516 … 143.9230469923935, 234.5189128148217, 1.0197882175752968, 0.19048212024807787, 1162.466061847129, 1168.8258601928706, 0.18666213351381167, 144453.56223654887, 10054.123067789495, -85.40932854551109]
[0.11608844570967189, 4.123269435043337e-5, 5.321779077597344e-5, 255.24937550869126, 287.19677065978607, 208.91246399469568, 31.407141823744052, 49.18129750484523, 243.7728973913202, 14.45064786272562 … 143.2148747590722, 233.48030821410865, 1.0176144773033036, 0.18763975074137146, 1151.0270265008958, 1156.7091554938063, 0.18401098194383467, 144456.3703720453, 10063.027094629098, -85.47042249594953]
[0.11615594839551835, 4.137267339973976e-5, 5.33988566967633e-5, 255.15822423386413, 287.14140165395537, 208.79069104335912, 31.444492103506605, 49.27554890457126, 243.63896944459316, 14.510313330128195 … 142.5266588778086, 232.4505787616062, 1.01528663154389, 0.18540560078986587, 1138.6926874652336, 1143.8358873415925, 0.18196508652491794, 144459.38988139568, 10072.117364448077, -85.51782755431448]
[0.11623831344571221, 4.151733614700252e-5, 5.3585987203993156e-5, 255.06472816214645, 287.08459944064526, 208.66599314853354, 31.482740294972167, 49.37249641868174, 243.50137910589646, 14.571794981919064 … 141.82430121812132, 231.38662596117442, 1.0126866807495154, 0.18354148259162642, 1125.2041482591612, 1129.8956684918414, 0.18029427904691805, 144462.7763236423, 10081.888795521523, -85.55715996221632]
[0.11633687605104438, 4.166507312070666e-5, 5.377710071790986e-5, 254.96998211671863, 287.02703071835845, 208.53988107046567, 31.52137114159804, 49.4709667161816, 243.36184478625816, 14.63435532218972 … 141.11394010942644, 230.29841067627248, 1.0098170633413128, 0.1819630483196959, 1110.9431576835855, 1115.2542442142521, 0.17891149019703545, 144466.53059020912, 10092.33729263554, -85.59049984206342]
[0.11634784946966706, 4.168020591875532e-5, 5.3796676985660515e-5, 254.9603195924157, 287.0211593505814, 208.52703558963728, 31.525301230914344, 49.48102035255478, 243.34761262828297, 14.640748694722598 … 141.04148135737876, 230.1866305267738, 1.0095103386271642, 0.18181516843656387, 1109.4645512207965, 1113.7402556728061, 0.17878353291084345, 144466.93289193546, 10093.437200922455, -85.59363341346629]
[0.11634784946966706, 4.168020591875532e-5, 5.3796676985660515e-5, 254.9603195924157, 287.0211593505814, 208.52703558963728, 31.525301230914344, 49.48102035255478, 243.34761262828297, 14.640748694722598 … 141.04148135737876, 230.1866305267738, 1.0095103386271642, 0.18181516843656387, 1109.4645512207965, 1113.7402556728061, 0.17878353291084345, 144466.93289193546, 10093.437200922455, -85.59363341346629]
@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.