diff --git a/GUI/backend/src/NTModelModule.jl b/GUI/backend/src/NTModelModule.jl index ce8b635..c68254f 100644 --- a/GUI/backend/src/NTModelModule.jl +++ b/GUI/backend/src/NTModelModule.jl @@ -19,7 +19,9 @@ mutable struct NTModel ps_units::Dict{Symbol, String} u0_desc::Dict{Symbol, String} ps_desc::Dict{Symbol, String} - + + # Parameters + # ---------- function NTModel() # Initial conditions @@ -27,13 +29,13 @@ mutable struct NTModel :GLN_e => 0.28, # mM :GLN_i => 0.07, # mM :GLU_e => 10.0, # mM - :GLU_i => 0.06, # mM - :GABA_i => 0.6, # mM - :GABA_v => 0.06, # mM + :GLU_i => 0.13, # mM + :GABA_i => 0.73, # mM + :GABA_v => 0.1, # mM :GLN_a => 3.8, # mM :GLU_a => 0.8, # mM :KET_b => 0.0, # mM - ]) + ]); # [X] = mM # Units for initial conditions u0_units = Dict([ @@ -46,7 +48,7 @@ mutable struct NTModel :GLN_a => "mM", :GLU_a => "mM", :KET_b => "mM", - ]) + ]); u0_desc = Dict([ :GLN_e => "Glutamine in excitatory neurons", @@ -58,29 +60,32 @@ mutable struct NTModel :GLN_a => "Glutamine in astrocytes", :GLU_a => "Glutamate in astrocytes", :KET_b => "Brain ketone concentration", - ]) + ]); # Parameters ps = Dict([ - # Neurotransmitter cycling + # Neurotranstmitter cycling :c_vpmas => 1.0, - :k_syn_e => 0.01912, - :k_syn_i => 0.79882, - :Vmax_GAD67 => 0.20769, - :KM_GAD67 => 0.2, - :Vmax_GT => 0.05860, - :KM_GT => 0.13359, - :Vmax_GAD65 => 0.05591, + :k_syn_e => 0.0191, + :k_syn_i => 0.4793, + :Vmax_GAD67 => 0.0959, + :KM_GAD67 => 0.13, + :Vmax_GT => 0.0582, + :KM_GT => 0.1563, + :Vmax_GAD65 => 0.0516, :KM_GAD65 => 0.01, + :Vmax_GDH => 0.096, + :KM_GDH => 0.8, + :h_GDH => 4.0, :Vmax_GS => 0.753, :KM_GS => 0.8, :KI_GS => 3.8, - :h_GS => 4.0, - :Vmax_GDH => 0.0764, - :KM_GDH => 0.8, - :h_GDH => 4.0, - :k_ox => 0.2, - :V_out => 0.012, + :V_anapl => 0.06, + :GLN_p => 0.5, + :GLN_ratio => 1/8, + :V_out0 => 0.012, + :Tmax_gln => 0.018, + :KT_gln => 0.25, # Energy metabolism :V_ATP => 11.94, # mM/min @@ -92,29 +97,32 @@ mutable struct NTModel :k_ket => 0.162, # 1/min :ATPket => 21.5, :KET_p => 0, - ]) + ]);# [V] = mM/min, [KM] = mM, [X] = mM # Units for parameters ps_units = Dict([ - # Neurotransmitter cycling + # Neurotranstmitter cycling :c_vpmas => "-", :k_syn_e => "1/min", - :k_syn_i => "1/min", + :k_syn_i => "1/min", :Vmax_GAD67 => "mM/min", :KM_GAD67 => "mM", :Vmax_GT => "mM/min", :KM_GT => "mM", :Vmax_GAD65 => "mM/min", :KM_GAD65 => "mM", - :Vmax_GS => "mM/min", - :KM_GS => "mM", - :KI_GS => "mM", - :h_GS => "-", :Vmax_GDH => "mM/min", :KM_GDH => "mM", :h_GDH => "-", - :k_ox => "-", - :V_out => "mM/min", + :Vmax_GS => "mM/min", + :KM_GS => "mM", + :KI_GS => "mM", + :V_anapl => "mM/min", + :GLN_p => "mM", + :GLN_ratio => "-", + :V_out0 => "mM/min", + :Tmax_gln => "mM/min", + :KT_gln => "mM", # Energy metabolism :V_ATP => "mM/min", @@ -126,28 +134,32 @@ mutable struct NTModel :k_ket => "1/min", :ATPket => "-", :KET_p => "mM", - ]) + + ]); ps_desc = Dict([ - # Neurotransmitter cycling + # Neurotranstmitter cycling :c_vpmas => "relative rate of PMAS", :k_syn_e => "rate constant for glutamate vesicle loading", - :k_syn_i => "rate constant for GABA vesicle loading", + :k_syn_i => "rate constant for GABA vesicle loading", :Vmax_GAD67 => "GAD67 maximum rate", :KM_GAD67 => "GAD67 Michaelis-Menten constant for glutamate", :Vmax_GT => "GT maximum rate", :KM_GT => "GT Michaelis-Menten constant for GABA", :Vmax_GAD65 => "GAD65 maximum rate", :KM_GAD65 => "GAD65 Michaelis-Menten constant for glutamate", + :Vmax_GDH => "GDH maximum rate", + :KM_GDH => "GDH Michaelis–Menten constant for glutamate", + :h_GDH => "GDH Hill exponent", :Vmax_GS => "GS maximum rate", :KM_GS => "GS Michaelis-Menten constant for glutamate", :KI_GS => "GS inhibition constant for glutamine", - :h_GS => "GS Hill exponent for the inhibition by GLN", - :Vmax_GDH => "GDH maximum rate", - :KM_GDH => "GDH Michaelis-Menten constant for glutamate", - :h_GDH => "GDH Hill exponent", - :k_ox => "relative portion of glutamate oxidation", - :V_out => "rate of GLN efflux", + :V_anapl => "rate of anaplerotic flux", + :GLN_p => "plasma glutamine concentration", + :GLN_ratio => "ratio of extracellular vs. intracellular glutamine", + :V_out0 => "baseline rate of glutamine efflux", + :Tmax_gln => "glutamine transport maximum rate", + :KT_gln => "glutamine transport Michaelis-Menten constant", # Energy metabolism :V_ATP => "rate of ATP synthesis", @@ -159,45 +171,44 @@ mutable struct NTModel :k_ket => "Ketone utilization rate constant", :ATPket => "ATP yield from ketones", :KET_p => "Plasma Ketone concentration", - ]) + ]); # Model definition + # ------------- + # Custom rate functions - mm_reversible(S1P, S1B, Vmax, Km1) = - Vmax / Km1 * (S1P - S1B) / (1 + S1P/Km1 + S1B/Km1) + mm_reversible(S1P, S1B, Vmax, Km) = Vmax * (S1P - S1B) / (Km + S1P + S1B); CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket) = ((V_ATP - ATPket*KET_b*k_ket)/ATPglc)/0.91 - Vpmas_e(cmrglc, c_vpmas) = c_vpmas*(cmrglc*0.91*0.75*0.83 - 0.077)/0.812 - Vpmas_i(cmrglc, c_vpmas) = c_vpmas*(cmrglc*0.91*0.75*0.17 - 0.011)/0.763 - + Vpmas_e(cmrglc) = (cmrglc*0.91*0.75*0.83 - 0.077)/0.812; + Vpmas_i(cmrglc) = (cmrglc*0.91*0.75*0.17 - 0.011)/0.763; + # Reaction network model = @reaction_network begin - - # Ketone transport + + # Ketone metabolism mm_reversible(KET_p, KET_b, Tmax_ket, KT_ket), 0 => KET_b, [description="Ketone transport"] - - # Ketone utilization k_ket*KET_b, KET_b => 0, [description="Ketone utilization"] # Glutamatergic neuron - Vpmas_e(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket), c_vpmas), GLN_e => GLU_e, [description="PMAS in e"] - k_syn_e * GLU_e, GLU_e => GLU_a, [description="GLU cycling"] + Vpmas_e(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket))*c_vpmas, GLN_e => GLU_e, [description="PMAS in e"] + k_syn_e * GLU_e, GLU_e => GLU_a, [description="GLU synaptic flux"] # GABAergic neuron - Vpmas_i(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket), c_vpmas), GLN_i => GLU_i, [description="PMAS in i"] + Vpmas_i(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket))*c_vpmas, GLN_i => GLU_i, [description="PMAS in i"] mm(GLU_i, Vmax_GAD67, KM_GAD67), GLU_i => GABA_i, [description="GAD67"] mm(GLU_i, Vmax_GAD65, KM_GAD65), GLU_i => GABA_v, [description="GAD65"] mm(GABA_i, Vmax_GT, KM_GT), GABA_i => GLU_i, [description="GT"] - k_syn_i * GABA_v, GABA_v => GLU_a, [description="GABA cycling"] + k_syn_i * GABA_v, GABA_v => GLU_a, [description="GABA synaptic flux"] # Astrocyte - k_ox * k_syn_e * GLU_e + V_out, 0 => GLU_a, [description="Anaplerosis"] + V_anapl, 0 => GLU_a, [description="Anaplerosis"] Vmax_GDH * (GLU_a/KM_GDH)^h_GDH / (1 + (GLU_a/KM_GDH)^h_GDH), GLU_a => 0, [description="GDH"] - Vmax_GS * GLU_a / (KM_GS*(1 + (GLN_a/KI_GS)^h_GS) + GLU_a), GLU_a => GLN_a, [description="GS"] - Vpmas_e(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket), c_vpmas), GLN_a => GLN_e, [description="GLN uptake in e"] - Vpmas_i(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket), c_vpmas), GLN_a => GLN_i, [description="GLN uptake in i"] - V_out, GLN_a => 0, [description="GLN efflux"] + Vmax_GS * GLU_a / (KM_GS*(1 + GLN_a/KI_GS) + GLU_a), GLU_a => GLN_a, [description="GS"] + Vpmas_e(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket))*c_vpmas, GLN_a => GLN_e, [description="GLN uptake in e"] + Vpmas_i(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket))*c_vpmas, GLN_a => GLN_i, [description="GLN uptake in i"] + V_out0 + mm_reversible(GLN_a*GLN_ratio, GLN_p, Tmax_gln, KT_gln), GLN_a => 0, [description="GLN efflux"] - end + end; return new(model, u0, ps, nothing, u0_units, ps_units, u0_desc, ps_desc) diff --git a/reproduce_paper/code/1_simulate_model 2.jl b/reproduce_paper/code/1_simulate_model 2.jl new file mode 100644 index 0000000..2c9a42e --- /dev/null +++ b/reproduce_paper/code/1_simulate_model 2.jl @@ -0,0 +1,438 @@ +""" + +@author: Botond B. Antal + +June 2026 + +This script runs the simulations and stores the results. + +""" + +using Catalyst +using OrdinaryDiffEq +using Printf +using DataFrames +using CSV +using Random +using HDF5 +using CairoMakie + +# Filepaths +OUTDIR = "" * "results/"; #TODO: insert your filepath here before results/ + +## +# Model +# ------------------------- + +# Burn-in time for simulations to ensure steady-state +burnin = 10000; + +# Custom rate functions +mm_reversible(S1P, S1B, Vmax, Km) = Vmax * (S1P - S1B) / (Km + S1P + S1B); +CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket) = ((V_ATP - ATPket*KET_b*k_ket)/ATPglc)/0.91 +Vpmas_e(cmrglc) = (cmrglc*0.91*0.75*0.83 - 0.077)/0.812; +Vpmas_i(cmrglc) = (cmrglc*0.91*0.75*0.17 - 0.011)/0.763; + +# Ketone input time-series +# KET_p(t) = 10.4*(exp(-(t-burnin)/89.2) - exp(-(t-burnin)/27.2))*1/(1+exp(-1e5* (t-burnin))) # For timecourse simulation +KET_p(t) = 0.0 # For MCA only + +# Create model +model = @reaction_network begin + + # Ketone metabolism + mm_reversible(KET_p(t), KET_b, Tmax_ket, KT_ket), 0 => KET_b, [description="Ketone transport"] + k_ket*KET_b, KET_b => 0, [description="Ketone utilization"] + + # Glutamatergic neuron + Vpmas_e(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket))*c_vpmas, GLN_e => GLU_e, [description="PMAS in e"] + k_syn_e * GLU_e, GLU_e => GLU_a, [description="GLU synaptic flux"] + + # GABAergic neuron + Vpmas_i(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket))*c_vpmas, GLN_i => GLU_i, [description="PMAS in i"] + mm(GLU_i, Vmax_GAD67, KM_GAD67), GLU_i => GABA_i, [description="GAD67"] + mm(GLU_i, Vmax_GAD65, KM_GAD65), GLU_i => GABA_v, [description="GAD65"] + mm(GABA_i, Vmax_GT, KM_GT), GABA_i => GLU_i, [description="GT"] + k_syn_i * GABA_v, GABA_v => GLU_a, [description="GABA synaptic flux"] + + # Astrocyte + V_anapl, 0 => GLU_a, [description="Anaplerosis"] + Vmax_GDH * (GLU_a/KM_GDH)^h_GDH / (1 + (GLU_a/KM_GDH)^h_GDH), GLU_a => 0, [description="GDH"] + Vmax_GS * GLU_a / (KM_GS*(1 + GLN_a/KI_GS) + GLU_a), GLU_a => GLN_a, [description="GS"] + Vpmas_e(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket))*c_vpmas, GLN_a => GLN_e, [description="GLN uptake in e"] + Vpmas_i(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket))*c_vpmas, GLN_a => GLN_i, [description="GLN uptake in i"] + V_out0 + mm_reversible(GLN_a*GLN_ratio, GLN_p, Tmax_gln, KT_gln), GLN_a => 0, [description="GLN efflux"] + + # Monitoring only (dummy reactions) + CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket), KET_b => KET_b, [description="CMRglc"] + CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket)*0.91*0.75, KET_b => KET_b, [description="CMRglc(ox)N"] + Vpmas_e(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket))*c_vpmas + + Vpmas_i(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket))*c_vpmas, + KET_b => KET_b, [description="VPMAS"] + +end; + +# Initial conditions +u0 = [ + :GLN_e => 0.28, # mM + :GLN_i => 0.07, # mM + :GLU_e => 10.0, # mM + :GLU_i => 0.13, # mM + :GABA_i => 0.73, # mM + :GABA_v => 0.1, # mM + :GLN_a => 3.8, # mM + :GLU_a => 0.8, # mM + :KET_b => 0.0, # mM + ]; # [X] = mM + +# Parameters +ps = Dict([ + + # Neurotranstmitter cycling + :k_syn_e => 0.0191, + :k_syn_i => 0.4793, + :Vmax_GAD67 => 0.0959, + :KM_GAD67 => 0.13, + :Vmax_GT => 0.0582, + :KM_GT => 0.1563, + :Vmax_GAD65 => 0.0516, + :KM_GAD65 => 0.01, + :Vmax_GDH => 0.096, + :KM_GDH => 0.8, + :h_GDH => 4.0, + :Vmax_GS => 0.753, + :KM_GS => 0.8, + :KI_GS => 3.8, + :V_anapl => 0.06, + :GLN_p => 0.5, + :GLN_ratio => 1/8, + :V_out0 => 0.012, + :Tmax_gln => 0.018, + :KT_gln => 0.25, + + # Energy metabolism + :V_ATP => 11.94, # mM/min + :ATPglc => 32.0, + + # Ketone metabolism + :Tmax_ket => 0.462, # mM/min + :KT_ket => 27.5, # mM + :k_ket => 0.162, # 1/min + :ATPket => 21.5, + + # MCA only + :c_vpmas => 1.0, # For MCA purposes only, must be set to 1. + + ]);# [V] = mM/min, [KM] = mM, [X] = mM + +# Time span +tmax = burnin + 500; +tspan = (0., tmax); + +# ODE problem +prob = ODEProblem(model, u0, tspan, ps); + +# Simulate +sol_raw = solve(prob, Vern7()); +sol = sol_raw[findall(sol_raw.t .>= burnin)]; # Cut off burn-in time + +## +# Single run of the model for testing +# ------------------------------------ + +# Make dataframe +df = DataFrame(sol); +rename!(df, [replace(item, "(t)" => "") for item in names(df)]); +df[!, :timestamp] = sol.t .- burnin; + +# Plot +inch=96; pt=4/3; +fig = Figure(size=(4.825inch, 3inch), fontsize=12pt) +ax = Axis(fig[1, 1]); +[lines!(ax, df[:, "timestamp"], df[:, label], label=label, linewidth=2) + for label in names(df)[2:end]]; + +CairoMakie.ylims!(0, 11.) +ax.xlabel = "Time"; +ax.ylabel = "Concentration"; +vlines!([60], color=:black, linestyle=:dash, linewidth=1); +Legend(fig[1, 2], ax, labelsize=8pt); +display(fig) + +# Print concentration changes +println("\nConcentrations:"); +for i in range(1, stop=length(unknowns(model))) + println(@sprintf("%s: %.4g -> %.4g", unknowns(model)[i], sol.u[1][i], sol.u[end][i])) +end; + +# Print reaction rates in steady state +println("\nReaction rates:"); +for i in range(1, stop=length(reactions(model))) + func_sym = reactionrates(model)[i]; + func = eval(build_function(func_sym, [p[1] for p in ps]..., :t, unknowns(model)...)); + rate = func([p[2] for p in ps]..., tmax, [sol[state][end] for state in unknowns(model)]...); + println(@sprintf("r=%.3g; ", rate) * Catalyst.getdescription(reactions(model)[i])) +end; + +# Compute fluxes for intermediate timepoints +items = []; +for i in range(1, stop=length(reactions(model))) + func_sym = reactionrates(model)[i]; + func = eval(build_function(func_sym, [p[1] for p in ps]..., :t, unknowns(model)...)); + reaction_name = Catalyst.getdescription(reactions(model)[i]); + rates = [func([p[2] for p in ps]..., tmax, [sol[state][j] for state in unknowns(model)]...) for j in 1:size(sol, 2)]; + push!(items, (Symbol(reaction_name) => rates)) +end; + +# Convert fluxes to dataframe +df_flux = DataFrame(items); +df_flux[!, :timestamp] = df.timestamp; + +# Save as CSV +CSV.write(OUTDIR * "model_timeseries.csv", df); +CSV.write(OUTDIR * "model_fluxes.csv", df_flux); + +## +# Full run with 1000 replicates +# ------------------------------------ + +# Set random seed for reproducibility +Random.seed!(1234); + +# Time span +tmax = burnin + 500; +tspan = (0., tmax); # [time] = min +save_t = range(0., stop=tmax, length=tmax+1); + +# Reference data +ref_start_glu = 10.93; +ref_start_gaba = 0.83; +ref_stdev_glu = 0.65 * ref_start_glu/6.89; +ref_stdev_gaba = 0.27 * ref_start_gaba/0.96; + +# Randomization settings +n_runs = 1000; # Number of replicates +frac = 0.10; # Fraction of variation for parameters (+/- %) + +# Open HDF5 file for writing +h5file = h5open(OUTDIR * "model_multirun.h5", "w"); + +# Try +try + + # Counter for successful runs + cnt_total = 0; + cnt_success = 0; + + # Create groups for pars, concs and fluxes + grp_info = create_group(h5file, "info") + grp_conc = create_group(h5file, "concentrations") + grp_flux = create_group(h5file, "fluxes") + grp_u0 = create_group(h5file, "u0") + grp_pars = create_group(h5file, "parameters") + + # Iterate through runs + while cnt_success < n_runs + + # Record attempt + cnt_total += 1; + + # Print status + print("\nRunning simulation #$cnt_total. Progress: ($(cnt_success)/$n_runs)") + + # Run simulation + # ---------- + + # Create new parameter set with random variations + ps_new = Dict( + k => v * (1 + frac * (2 * rand() - 1)) for (k, v) in ps); + + # Create a new set of initial concentrations + u0_new = Dict( + k => v * (1 + frac * (2 * rand() - 1)) for (k, v) in u0); + + # Fix parameters related to activation state (out-of-scope for model) + ps_new[:c_vpmas] = ps[:c_vpmas] + + # Remake problem + prob = remake(prob; u0=u0_new, p=ps_new, tspan=tspan); + + # Simulate + sol_raw = solve(prob, Vern7(), saveat=save_t); + + # Check simulation results + # ------------------ + + # Check for steady state + mask = abs.(sol_raw(burnin-10, Val{1})) .> 1e-4; + if sum(mask) > 0 + print(" // Not in steady-state: ($(sum(mask))) ", + join(unknowns(model)[mask], ", ")); + continue + end + + # Get starting state once steady-state is reached + start_state = sol_raw(burnin); + start_glu = start_state[3] + start_state[4] + start_state[6] + start_gaba = start_state[7] + start_state[8] + + # Check if starting state is within reference range + mask = [ + abs(start_glu - ref_start_glu) > ref_stdev_glu, + abs(start_gaba - ref_start_gaba) > ref_stdev_gaba + ] + if sum(mask) > 0 + print(" // Starting state outside reference range: ", + join(["GLU", "GABA"][mask], ", ")); + continue + end + + # Cut off burn-in time + sol = sol_raw[findall(save_t .>= burnin)]; + + # Extract and refine varnames + var_names = replace.(string.(unknowns(model)), "(t)" => "") + + # Parameters + # --------------- + grp = create_group(grp_pars, "run_$(cnt_success+1)") + grp["par_names"] = [string(k) for k in keys(ps_new)] + grp["pars"] = [v for v in values(ps_new)] + + # Concentrations + # -------------- + + # Save + grp = create_group(grp_conc, "run_$(cnt_success+1)") + grp["var_names"] = var_names + grp["t"] = sol.t .- burnin + grp["concs"] = Array(sol) + + # u0 + # ---------- + + # Save + grp = create_group(grp_u0, "run_$(cnt_success+1)") + grp["var_names"] = var_names + grp["u0_names"] = [string(k) for k in keys(u0_new)] + grp["u0"] = [v for v in values(u0_new)] + + # Fluxes + # ------------- + + # Initialize collection for columns of rates + items = []; + + # Iterate through reactions and extract the reaction description and rates + for i in range(1, stop=length(reactions(model))) + func_sym = reactionrates(model)[i]; + func = eval(build_function(func_sym, [p[1] for p in ps]..., :t, unknowns(model)...)); + reaction_name = Catalyst.getdescription(reactions(model)[i]); + rates = [func([p[2] for p in ps]..., tmax, [sol[state][j] for state in unknowns(model)]...) for j in 1:size(sol, 2)]; + push!(items, (Symbol(reaction_name) => rates)) + end; + + # Save + grp = create_group(grp_flux, "run_$(cnt_success+1)") + grp["reaction_names"] = [string(item[1]) for item in items] + grp["t"] = sol.t .- burnin + grp["rates"] = hcat([item[2] for item in items]...) # Convert to matrix + + # Mark run as successful + print(" // Successful run.") + cnt_success += 1; + + end + + # Count successful runs + println("\nSuccessful runs: $cnt_success out of $cnt_total.") + +finally + + # Close h5 files + close(h5file) + +end + +## +# MCA - concentration control +# ------------------------- + +""" +Rerun the model definition with "KET_p(t) = 0.0" before running this section! +""" + +# Setup +C_rows = [] # Collection for control coefficients +ϵ = 0.01 # Extent of perturbation +tspan = (0., burnin); # Longer time span to reach steady state + +# Iterate through the parameters +for par in keys(ps) + + # Reset the parameter set + ps_test = copy(ps) + + # Center + ps_test[par] = ps[par]; + prob = remake(prob; u0=u0, p=ps_test, tspan=tspan); + sol = solve(prob, Vern7()); + S0 = [sol[state][end] for state in unknowns(model)] + + # Increment 1 + ps_test[par] = ps[par] * (1 - ϵ); + prob = remake(prob; u0=u0, p=ps_test, tspan=tspan); + sol = solve(prob, Vern7()); + S1 = [sol[state][end] for state in unknowns(model)] + + # Increment 2 + ps_test[par] = ps[par] * (1 + ϵ); + prob = remake(prob; u0=u0, p=ps_test, tspan=tspan); + sol = solve(prob, Vern7()); + S2 = [sol[state][end] for state in unknowns(model)] + + # Compute the control coefficient based on central differences + C = (S2 .- S1) ./ (S0 .* 2ϵ) + push!(C_rows, C) + +end + +# Compile results into dataframe +df = DataFrame(C_rows, :auto) +rename!(df, String.(keys(ps))) + +# Add state names +df[!, :state] = string.(unknowns(model)); + +# Save as CSV +CSV.write(OUTDIR * "model_mca_concentration.csv", df); + +# Create heatmap +mat = Matrix(df[!, 1:end-1])' +fig = Figure(size=(5.5inch, 4.5inch), fontsize=12pt) #, backgroundcolor=:transparent); +ax = Axis(fig[1, 1], + aspect=DataAspect(), + title="Concentration control coefficients", + xticks = (1:size(df, 2)-1, names(df)[1:end-1]), + yticks = (1:size(df, 1), df[!, :state]), + xticklabelrotation = pi / 2., + xgridvisible = false, + ygridvisible = false) +hm = heatmap!(ax, mat, colormap=:seismic, colorrange =(-30, +30),) +Colorbar(fig[1, 2], hm) + +# Add text labels to the heatmap +for i in 1:size(mat, 1), j in 1:size(mat, 2) + text!(ax, i, j, text=@sprintf("%.1f", mat[i, j]), + align=(:center, :center), color=:black, fontsize=8pt) +end + +# Add grid lines manually +for x in 1:size(mat, 1) + lines!(ax, [x-0.5, x-0.5], [0, size(mat, 2)+0.5], color=:black, linewidth=0.5) +end +for y in 1:size(mat, 2) + lines!(ax, [0, size(mat, 1)+0.5], [y-0.5, y-0.5], color=:black, linewidth=0.5) +end + +display(fig) diff --git a/reproduce_paper/code/1_simulate_model.jl b/reproduce_paper/code/1_simulate_model.jl deleted file mode 100644 index 2fabd97..0000000 --- a/reproduce_paper/code/1_simulate_model.jl +++ /dev/null @@ -1,262 +0,0 @@ -""" - -@author: Botond B. Antal - -January 2026 - -This script runs the simulations and saves the results. - -""" - -using Catalyst -using OrdinaryDiffEq -using Printf -using DataFrames -using CSV -using CairoMakie - -OUTDIR = "" * "results/"; #TODO: insert your filepath here before results/ - -## -# Model -# ------------------------- - -# Custom rate functions -mm_reversible(S1P, S1B, Vmax, Km1) = - Vmax / Km1 * (S1P - S1B) / (1 + S1P/Km1 + S1B/Km1); -CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket) = ((V_ATP - ATPket*KET_b*k_ket)/ATPglc)/0.91 -Vpmas_e(cmrglc, c_vpmas) = c_vpmas*(cmrglc*0.91*0.75*0.83 - 0.077)/0.812; -Vpmas_i(cmrglc, c_vpmas) = c_vpmas*(cmrglc*0.91*0.75*0.17 - 0.011)/0.763; - -# Ketone input time-series -KET_p(t) = 10.4*(exp(-(t)/89.2) - exp(-t/27.2)) # For comparison with test dataset -# KET_p(t) = 0.0 # For MCA - -# Create model -model = @reaction_network begin - - # Ketone transport - mm_reversible(KET_p(t), KET_b, Tmax_ket, KT_ket), 0 => KET_b, [description="Ketone transport"] - - # Ketone utilization - k_ket*KET_b, KET_b => 0, [description="Ketone utilization"] - - # Glutamatergic neuron - Vpmas_e(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket), c_vpmas), GLN_e => GLU_e, [description="PMAS in e"] - k_syn_e * GLU_e, GLU_e => GLU_a, [description="GLU cycling"] - - # GABAergic neuron - Vpmas_i(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket), c_vpmas), GLN_i => GLU_i, [description="PMAS in i"] - mm(GLU_i, Vmax_GAD67, KM_GAD67), GLU_i => GABA_i, [description="GAD67"] - mm(GLU_i, Vmax_GAD65, KM_GAD65), GLU_i => GABA_v, [description="GAD65"] - mm(GABA_i, Vmax_GT, KM_GT), GABA_i => GLU_i, [description="GT"] - k_syn_i * GABA_v, GABA_v => GLU_a, [description="GABA cycling"] - - # Astrocyte - k_ox * k_syn_e * GLU_e + V_out, 0 => GLU_a, [description="Anaplerosis"] - Vmax_GDH * (GLU_a/KM_GDH)^h_GDH / (1 + (GLU_a/KM_GDH)^h_GDH), GLU_a => 0, [description="GDH"] - Vmax_GS * GLU_a / (KM_GS*(1 + (GLN_a/KI_GS)^h_GS) + GLU_a), GLU_a => GLN_a, [description="GS"] - Vpmas_e(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket), c_vpmas), GLN_a => GLN_e, [description="GLN uptake in e"] - Vpmas_i(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket), c_vpmas), GLN_a => GLN_i, [description="GLN uptake in i"] - V_out, GLN_a => 0, [description="GLN efflux"] - - # Monitoring - CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket), KET_b => KET_b, [description="CMRglc"] - CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket)*0.91*0.75, KET_b => KET_b, [description="CMRglc(ox)N"] - Vpmas_e(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket), c_vpmas) + - Vpmas_i(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket), c_vpmas), - KET_b => KET_b, [description="VPMAS"] - -end; - -# Initial conditions -u0 = [ - :GLN_e => 0.28, # mM - :GLN_i => 0.07, # mM - :GLU_e => 10.0, # mM - :GLU_i => 0.06, # mM - :GABA_i => 0.6, # mM - :GABA_v => 0.06, # mM - :GLN_a => 3.8, # mM - :GLU_a => 0.8, # mM - :KET_b => 0.0, # mM - ]; # [X] = mM - -# Parameters -ps = Dict([ - - # Neurotranstmitter cycling - :c_vpmas => 1.0, - :k_syn_e => 0.01912, - :k_syn_i => 0.79882, - :Vmax_GAD67 => 0.20769, - :KM_GAD67 => 0.2, - :Vmax_GT => 0.05860, - :KM_GT => 0.13359, - :Vmax_GAD65 => 0.05591, - :KM_GAD65 => 0.01, - :Vmax_GS => 0.753, - :KM_GS => 0.8, - :KI_GS => 3.8, - :h_GS => 4.0, - :Vmax_GDH => 0.0764, - :KM_GDH => 0.8, - :h_GDH => 4.0, - :k_ox => 0.2, - :V_out => 0.012, - - # Energy metabolism - :V_ATP => 11.94, # mM/min - :ATPglc => 32.0, - - # Ketone metabolism - :Tmax_ket => 0.462, # mM/min - :KT_ket => 27.5, # mM - :k_ket => 0.162, # 1/min - :ATPket => 21.5, - - ]);# [V] = mM/min, [KM] = mM, [X] = mM - -# Time span -tmax = 1000; -tspan = (0., tmax); - -# ODE problem -prob = ODEProblem(model, u0, tspan, ps); - -# Simulate -sol = solve(prob, Vern7()); - -# Analyze simulation results -# ---------- - -# Make dataframe -df = DataFrame(sol); -rename!(df, [replace(item, "(t)" => "") for item in names(df)]); - -# Plot -inch=96; pt=4/3; -fig = Figure(size=(4.825inch, 3inch), fontsize=12pt) -ax = Axis(fig[1, 1]); -[lines!(ax, df[:, "timestamp"], df[:, label], label=label, linewidth=2) - for label in names(df)[2:end]]; - -CairoMakie.ylims!(0, 11.) -ax.xlabel = "Time"; -ax.ylabel = "Concentration"; -vlines!([60], color=:black, linestyle=:dash, linewidth=1); -Legend(fig[1, 2], ax, labelsize=8pt); -display(fig) - -# Print concentration changes -println("\nConcentrations:"); -for i in range(1, stop=length(unknowns(model))) - println(@sprintf("%s: %.4g -> %.4g", unknowns(model)[i], sol.u[1][i], sol.u[end][i])) -end; - -# Print reaction rates in steady state -println("\nReaction rates:"); -for i in range(1, stop=length(reactions(model))) - func_sym = reactionrates(model)[i]; - func = eval(build_function(func_sym, [p[1] for p in ps]..., :t, unknowns(model)...)); - rate = func([p[2] for p in ps]..., tmax, [sol[state][end] for state in unknowns(model)]...); - println(@sprintf("r=%.3g; ", rate) * Catalyst.getdescription(reactions(model)[i])) -end; - - -# Compute fluxes for intermediate timepoints -items = []; -for i in range(1, stop=length(reactions(model))) - func_sym = reactionrates(model)[i]; - func = eval(build_function(func_sym, [p[1] for p in ps]..., :t, unknowns(model)...)); - reaction_name = Catalyst.getdescription(reactions(model)[i]); - rates = [func([p[2] for p in ps]..., tmax, [sol[state][j] for state in unknowns(model)]...) for j in 1:size(sol, 2)]; - push!(items, (Symbol(reaction_name) => rates)) -end; - -# Convert fluxes to dataframe -df_flux = DataFrame(items); -df_flux[!, :timestamp] = df.timestamp; - -# Save as CSV -CSV.write(OUTDIR * "model_timeseries.csv", df); -CSV.write(OUTDIR * "model_fluxes.csv", df_flux); - -## -# MCA - concentration control -# ------------------------- - -# Setup -C_rows = [] # Collection for control coefficients -ϵ = 0.01 # Extent of perturbation -tspan = (0., 1000.); # Longer time span to reach steady state - -# Iterate through the parameters -for par in keys(ps) - - # Reset the parameter set - ps_test = copy(ps) - - # Center - ps_test[par] = ps[par]; - prob = remake(prob; u0=u0, p=ps_test, tspan=tspan); - sol = solve(prob, Vern7()); - S0 = [sol[state][end] for state in unknowns(model)] - - # Increment 1 - ps_test[par] = ps[par] * (1 - ϵ); - prob = remake(prob; u0=u0, p=ps_test, tspan=tspan); - sol = solve(prob, Vern7()); - S1 = [sol[state][end] for state in unknowns(model)] - - # Increment 2 - ps_test[par] = ps[par] * (1 + ϵ); - prob = remake(prob; u0=u0, p=ps_test, tspan=tspan); - sol = solve(prob, Vern7()); - S2 = [sol[state][end] for state in unknowns(model)] - - # Compute the control coefficient based on central differences - C = (S2 .- S1) ./ (S0 .* 2ϵ) - push!(C_rows, C) - -end - -# Compile results into dataframe -df = DataFrame(C_rows, :auto) -rename!(df, String.(keys(ps))) - -# Add state names -df[!, :state] = string.(unknowns(model)); - -# Save as CSV -CSV.write(OUTDIR * "model_mca_concentration.csv", df); - -# Create heatmap -mat = Matrix(df[!, 1:end-1])' -fig = Figure(size=(5.5inch, 4.5inch), fontsize=12pt) #, backgroundcolor=:transparent); -ax = Axis(fig[1, 1], - aspect=DataAspect(), - title="Concentration control coefficients", - xticks = (1:size(df, 2)-1, names(df)[1:end-1]), - yticks = (1:size(df, 1), df[!, :state]), - xticklabelrotation = pi / 2., - xgridvisible = false, - ygridvisible = false) -hm = heatmap!(ax, mat, colormap=:seismic, colorrange =(-600, +600),) -Colorbar(fig[1, 2], hm) - -# Add text labels to the heatmap -for i in 1:size(mat, 1), j in 1:size(mat, 2) - text!(ax, i, j, text=@sprintf("%.1f", mat[i, j]), - align=(:center, :center), color=:black, fontsize=8pt) -end - -# Add grid lines manually -for x in 1:size(mat, 1) - lines!(ax, [x-0.5, x-0.5], [0, size(mat, 2)+0.5], color=:black, linewidth=0.5) -end -for y in 1:size(mat, 2) - lines!(ax, [0, size(mat, 1)+0.5], [y-0.5, y-0.5], color=:black, linewidth=0.5) -end - -display(fig) diff --git a/reproduce_paper/code/2_plot_results 2.py b/reproduce_paper/code/2_plot_results 2.py new file mode 100644 index 0000000..6936599 --- /dev/null +++ b/reproduce_paper/code/2_plot_results 2.py @@ -0,0 +1,784 @@ +""" + +@author: Botond B. Antal + +June 2026 + +This script plots the simulated results. + +""" + +# %% + +import os +import pandas as pd +import numpy as np +import matplotlib.pyplot as plt +import seaborn as sns +import h5py + +# %% +# =================================================================== +# Setup +# =================================================================== + +# Directories +OUTDIR = "" + "results/" #TODO: insert your filepath here before results/ + +# rcParams +plt.rcParams.update({'font.size': 10}) + +# Experiment data +t_meas = 67 +t_meas_error = 4.1 +exp_means = {"GLU": -10.6, "GABA": -36.9} +exp_errors = {"GLU": 1.2*1.96, "GABA": 4.2*1.96} # 95% CI! +exp_scaling = {"GLU": 11, "GABA": 0.83} + +# %% +# =================================================================== +# Load simulation results with replicates +# =================================================================== + +# Initialize dataframe +df_list_params = [] +df_list_conc = [] +df_list_conc_norm = [] +df_list_flux = [] +df_list_flux_norm = [] +df_list_u0 = [] + +# Open conc data and make long format dataframe for each run +with h5py.File(OUTDIR + "model_multirun.h5", "r") as f: + + # Concentration group + grp_conc = f["concentrations"] + + # Get run names + run_names = list(grp_conc.keys()) + + # Iterate over runs + for run in run_names: + + # Extract data + t = grp_conc[run]["t"][:] + var_names = [item.decode() for item in grp_conc[run]["var_names"]] + Y = grp_conc[run]["concs"][:] + + # Create dataframe for this run + df_run = pd.DataFrame(Y, columns=var_names) + df_run["timestamp"] = t + df_run["run"] = run + + # Append to list + df_list_conc.append(df_run) + + # Normalize with first timepoint and convert to percent change + df_run_norm = df_run.copy() + for var in var_names: + df_run_norm[var] = 100 * (df_run[var] - df_run[var][0]) / df_run[var][0] + df_run_norm["timestamp"] = t + + # Add combined columns + # Glutamate + cols = ["GLU_e", "GLU_i", "GLU_a"] + df_run_norm["GLU_tot"] = 100 * (df_run[cols].sum(axis=1) - df_run[cols].sum(axis=1)[0]) / \ + df_run[cols].sum(axis=1)[0] + + # GABA + cols = ["GABA_i", "GABA_v"] + df_run_norm["GABA_tot"] = 100 * (df_run[cols].sum(axis=1) - df_run[cols].sum(axis=1)[0]) / \ + df_run[cols].sum(axis=1)[0] + + # Append normalized dataframe to list + df_list_conc_norm.append(df_run_norm) + + # Flux group + grp_flux = f["fluxes"] + + # Run names + run_names = list(grp_flux.keys()) + + for run in run_names: + + # Extract data + t = grp_flux[run]["t"][:] + reaction_names = [item.decode() for item in grp_flux[run]["reaction_names"]] + Y = grp_flux[run]["rates"][:] # (reactions, time) + + # Create dataframe for this run + df_run = pd.DataFrame(Y.T, columns=reaction_names) + df_run["timestamp"] = t + df_run["run"] = run + + # Append to list + df_list_flux.append(df_run) + + # Normalize with first timepoint and convert to percent change + df_run_norm = df_run.copy() + for reaction in reaction_names: + df_run_norm[reaction] = 100 * (df_run[reaction] - df_run[reaction][0]) / df_run[reaction][0] + df_run_norm["timestamp"] = t + + # Append normalized dataframe to list + df_list_flux_norm.append(df_run_norm) + + # Parameter group + grp_params = f["parameters"] + + # Run names + run_names = list(grp_params.keys()) + + # Iterate over runs + for run in run_names: + + # Extract data + par_names = [item.decode() for item in grp_params[run]["par_names"]] + values = grp_params[run]["pars"][:] + + # Create dataframe for this run + df_run = pd.DataFrame(values[:, None].T, columns=par_names) + df_run["run"] = run + + # Append to list + df_list_params.append(df_run) + + # Initial conditions group + grp_u0 = f["u0"] + + # Run names + run_names = list(grp_u0.keys()) + + # Iterate over runs + for run in run_names: + + # Extract data + var_names = [item.decode() for item in grp_u0[run]["var_names"]] + values = grp_u0[run]["u0"][:] + + # Create dataframe for this run + df_run = pd.DataFrame(values[:, None].T, columns=var_names) + df_run["run"] = run + + # Append to list + df_list_u0.append(df_run) + +# Concat dataframes +df_params = pd.concat(df_list_params, ignore_index=True) +df_conc = pd.concat(df_list_conc, ignore_index=True) +df_conc_norm = pd.concat(df_list_conc_norm, ignore_index=True) +df_flux = pd.concat(df_list_flux, ignore_index=True) +df_flux_norm = pd.concat(df_list_flux_norm, ignore_index=True) +df_u0 = pd.concat(df_list_u0, ignore_index=True) + +# %% +# =================================================================== +# Describe the results +# =================================================================== + +# Take the mean across runs for each timepoint +df_flux_mean = df_flux.drop("run", axis=1).groupby("timestamp").mean().reset_index() +df_conc_mean = df_conc.drop("run", axis=1).groupby("timestamp").mean().reset_index() +df_conc_norm_mean = df_conc_norm.drop("run", axis=1).groupby("timestamp").mean().reset_index() + +# Relative concentrations at t_meas +y_tmeas = df_conc_norm_mean.iloc[df_conc_norm_mean['timestamp'].sub(t_meas).abs().idxmin()] + +# Print values +print(f"CMRglc starts from {df_flux_mean['CMRglc'][0]:.3f} mM/min " + f"and is lowest at {df_flux_mean['CMRglc'].min():.3f} mM/min " + f"which occurs at t={df_flux_mean['timestamp'][df_flux_mean['CMRglc'].idxmin()]:.2f}) min") +print(f"CMRket starts from {df_flux_mean['Ketone utilization'][0]:.3f} mM/min " + f"and is highest at {df_flux_mean['Ketone utilization'].max():.3f} mM/min " + f"which occurs at t={df_flux_mean['timestamp'][df_flux_mean['Ketone utilization'].idxmax()]:.2f} min") +print(f"PMAS_e starts from {df_flux_mean['PMAS in e'][0]:.3f} mM/min " + f"and is lowest at {df_flux_mean['PMAS in e'].min():.3f} mM/min " + f"which occurs at t={df_flux_mean['timestamp'][df_flux_mean['PMAS in e'].idxmin()]:.2f} min") +print(f"PMAS_i starts from {df_flux_mean['PMAS in i'][0]:.3f} mM/min " + f"and is lowest at {df_flux_mean['PMAS in i'].min():.3f} mM/min " + f"which occurs at t={df_flux_mean['timestamp'][df_flux_mean['PMAS in i'].idxmin()]:.2f} min") +print(f"VGDH starts from {df_flux_mean['GDH'][0]:.3f} mM/min " + f"and is lowest at {df_flux_mean['GDH'].min():.3f} mM/min " + f"which occurs at t={df_flux_mean['timestamp'][df_flux_mean['GDH'].idxmin()]:.2f} min") +print(f"VGS starts from {df_flux_mean['GS'][0]:.3f} mM/min " + f"and is lowest at {df_flux_mean['GS'].min():.3f} mM/min " + f"which occurs at t={df_flux_mean['timestamp'][df_flux_mean['GS'].idxmin()]:.2f} min") +print(f"VGAD65 starts from {df_flux_mean['GAD65'][0]:.3f} mM/min " + f"and is lowest at {df_flux_mean['GAD65'].min():.3f} mM/min " + f"which occurs at t={df_flux_mean['timestamp'][df_flux_mean['GAD65'].idxmin()]:.2f} min") +print(f"VGAD67 starts from {df_flux_mean['GAD67'][0]:.3f} mM/min " + f"and is lowest at {df_flux_mean['GAD67'].min():.3f} mM/min " + f"which occurs at t={df_flux_mean['timestamp'][df_flux_mean['GAD67'].idxmin()]:.2f} min") +print(f"VGT starts from {df_flux_mean['GT'][0]:.3f} mM/min " + f"and is lowest at {df_flux_mean['GT'].min():.3f} mM/min " + f"which occurs at t={df_flux_mean['timestamp'][df_flux_mean['GT'].idxmin()]:.2f} min") + +print("\n####\n") +# Print values +print(f"GLU_e starts from {df_conc_mean['GLU_e'][0]:.3f} mM " + f"and is lowest at {df_conc_mean['GLU_e'].min():.3f} mM " + f"(-{100 * (df_conc_mean['GLU_e'][0] - df_conc_mean['GLU_e'].min()) / df_conc_mean['GLU_e'][0]:.1f}%) " + f"which occurs at t={df_conc_mean['timestamp'][df_conc_mean['GLU_e'].idxmin()]:.2f} min") +print(f"GLU_a starts from {df_conc_mean['GLU_a'][0]:.3f} mM " + f"and is lowest at {df_conc_mean['GLU_a'].min():.3f} mM " + f"(-{100 * (df_conc_mean['GLU_a'][0] - df_conc_mean['GLU_a'].min()) / df_conc_mean['GLU_a'][0]:.1f}%) " + f"which occurs at t={df_conc_mean['timestamp'][df_conc_mean['GLU_a'].idxmin()]:.2f} min") +print(f"GLU_a starts from {df_conc_mean['GLU_a'][0]:.3f} mM " + f"and is highest at {df_conc_mean['GLU_a'].max():.3f} mM " + f"(-{100 * (df_conc_mean['GLU_a'][0] - df_conc_mean['GLU_a'].max()) / df_conc_mean['GLU_a'][0]:.1f}%) " + f"which occurs at t={df_conc_mean['timestamp'][df_conc_mean['GLU_a'].idxmax()]:.2f} min") +print(f"GLN_a starts from {df_conc_mean['GLN_a'][0]:.3f} mM " + f"and is highest at {df_conc_mean['GLN_a'].max():.3f} mM " + f"(-{100 * (df_conc_mean['GLN_a'][0] - df_conc_mean['GLN_a'].max()) / df_conc_mean['GLN_a'][0]:.1f}%) " + f"which occurs at t={df_conc_mean['timestamp'][df_conc_mean['GLN_a'].idxmax()]:.2f} min") +print(f"GABA_i starts from {df_conc_mean['GABA_i'][0]:.3f} mM " + f"and is lowest at {df_conc_mean['GABA_i'].min():.3f} mM " + f"(-{100 * (df_conc_mean['GABA_i'][0] - df_conc_mean['GABA_i'].min()) / df_conc_mean['GABA_i'][0]:.1f}%) " + f"which occurs at t={df_conc_mean['timestamp'][df_conc_mean['GABA_i'].idxmin()]:.2f} min") +print(f"GABA_v starts from {df_conc_mean['GABA_v'][0]:.3f} mM " + f"and is lowest at {df_conc_mean['GABA_v'].min():.3f} mM " + f"(-{100 * (df_conc_mean['GABA_v'][0] - df_conc_mean['GABA_v'].min()) / df_conc_mean['GABA_v'][0]:.1f}%) " + f"which occurs at t={df_conc_mean['timestamp'][df_conc_mean['GABA_v'].idxmin()]:.2f} min") +print(f"KET_b starts from {df_conc_mean['KET_b'][0]:.3f} mM " + f"and is highest at {df_conc_mean['KET_b'].max():.3f} mM " + f"(-{100 * (df_conc_mean['KET_b'][0] - df_conc_mean['KET_b'].max()) / df_conc_mean['KET_b'][0]:.1f}%) " + f"which occurs at t={df_conc_mean['timestamp'][df_conc_mean['KET_b'].idxmax()]:.2f} min") + +# Print total (relative) changes at t_meas +print(f"Glutamate: {y_tmeas['GLU_tot']:.1f}%") +print(f"GABA: {y_tmeas['GABA_tot']:.1f}%") + +# %% +# =================================================================== +# Figure 2: ketone metabolism +# =================================================================== + +# Load data +data_bloodbhb_datapoints = pd.read_csv(OUTDIR + "martin-arrowsmith_etal_blood_bhb_timecourse.csv") +data_bloodbhb_timeseries = pd.read_csv(OUTDIR + f"model_blood_bhb_input_timeseries.csv") +data_brainbhb_datapoints = pd.read_csv(OUTDIR + "mujica-parodi_etal_brain_bhb_timeseries.csv") + +# Take the max of the mean ketone time-course across runs +ket_b_max = df_conc.groupby("timestamp")["KET_b"].mean().max() + +# Figure +plt.figure(figsize=(8.25, 3.3)) + +# Subplot 1 +plt.subplot(1, 3, 1) +plt.scatter( + data_bloodbhb_datapoints['time'], + data_bloodbhb_datapoints['concentration'], + label='Measured\n(Martin-\nArrowsmith\net al.)', + color='black', s=80, zorder=2) +plt.plot( + data_bloodbhb_timeseries['time'], + data_bloodbhb_timeseries['concentration'], + lw=1.5, label='Fitted (double\nexponential)', + color='crimson', zorder=3) +plt.xlabel('Time after bolus (min)') +plt.ylabel(r'Blood D-$\beta$HB concentration (mM)') +# plt.title(r"Blood D-$\beta$HB") +plt.axhline(0, color='gray', lw=1, ls='-') +plt.xlim(0, 220) +plt.xticks([0, 50, 100, 150, 200]) +leg = plt.legend(framealpha=1, fontsize=8, bbox_to_anchor=(0.08, 0.2), + edgecolor="black", loc="center left") +plt.grid(alpha=0.4, zorder=1) +leg.set_in_layout(False) + +# Subplot 2 +plt.subplot(1, 3, 2) +plt.scatter( + data_brainbhb_datapoints['x'], + data_brainbhb_datapoints['y']/12, + label='Measured and rescaled\n(Mujica-Parodi et al.)', + color='black', s=80, zorder=2) +sns.lineplot(data=df_conc, x="timestamp", y="KET_b", + errorbar=('sd'), label='Simulated', lw=1.5, color='dodgerblue') +plt.xlabel('Time after bolus (min)') +plt.ylabel(r'Brain D-$\beta$HB concentration (mM)') +# plt.title(r"Brain D-$\beta$HB") +plt.axhline(0, color='gray', lw=1, ls='-') +plt.xlim(0, 220) +plt.xticks([0, 50, 100, 150, 200]) +plt.legend(framealpha=1, fontsize=8, edgecolor="black") +plt.grid(alpha=0.4, zorder=1) + +# Subplot 3 +plt.subplot(1, 3, 3) +plt.plot( + data_bloodbhb_timeseries['time'], + data_bloodbhb_timeseries['concentration']/data_bloodbhb_timeseries['concentration'].max(), + lw=1.5, label='Blood (input)', color='crimson', zorder=2) +sns.lineplot(x=df_conc["timestamp"], y=df_conc["KET_b"]/ket_b_max, + errorbar=('sd'), label='Brain (simulated)', lw=1.5, color='dodgerblue') +plt.xlabel('Time after bolus (min)') +plt.ylabel('Concentration (relative)') +# plt.title("Normalized\n" + r'Blood & Brain D-$\beta$HB') +plt.axhline(0, color='gray', lw=1, ls='-') +plt.xlim([0, 100]) +plt.ylim([plt.gca().get_ylim()[0], 1.12]) +plt.legend(framealpha=1, fontsize=8, edgecolor="black") +plt.grid(alpha=0.4, zorder=0) + +# Save +plt.tight_layout() +plt.savefig(OUTDIR + f"figs/fig_ketone_pk.pdf", dpi=300, transparent=True) + +# %% +# =================================================================== +# Figure 3: neurotransmitter metabolism +# ================================================================== + +# Figure +plt.figure(figsize=(8.25, 3.3)) + +# Fluxes (A) +# -------------- + +# Subplot +plt.subplot(1, 4, 1) + +# Plot +sns.lineplot(data=df_flux, x="timestamp", y="CMRglc", errorbar=('sd'), + label='CMR$_{glc}$', lw=1.5, color="#000000") +sns.lineplot(data=df_flux, x="timestamp", y="Ketone utilization", errorbar=('sd'), + label='CMR$_{Ket}$', lw=1.5, color="mediumvioletred") +sns.lineplot(data=df_flux, x="timestamp", y="GS", errorbar=('sd'), + label='V$_{GS}$', lw=1.5, linestyle='-', color="#44AB42F7") +sns.lineplot(data=df_flux, x="timestamp", y="PMAS in e", errorbar=('sd'), + label='V$_{PMAS}^e$', lw=1.5, color='#0066FF') +sns.lineplot(data=df_flux, x="timestamp", y="PMAS in i", errorbar=('sd'), + label='V$_{PMAS}^i$', lw=1.5, color="#D60000F7") + + +# Format +plt.axhline(0, color='gray', lw=1, ls='-') +plt.xlabel('Time after bolus (min)') +plt.ylabel('Flux (mM/min)', labelpad=1) +plt.xticks([0, 50, 100, 150, 200]) +plt.xlim(0, 220) +plt.ylim([-0.02, 0.55]) +plt.grid(alpha=0.4) +# plt.title("Fluxes", fontsize=10) +leg = plt.legend(framealpha=1, fontsize=7, edgecolor="black", + bbox_to_anchor=(0.08, 0.86), loc="center left", + ncol=2, columnspacing=0.3, handlelength=1) +leg.set_in_layout(False) + +# Compartmental concentrations (B) +# -------------- + +# Subplot +plt.subplot(1, 4, 2) + +# Plot +sns.lineplot(data=df_conc, x="timestamp", y="GLU_e", + errorbar=('sd'), label='[Glu]$_{e}$', lw=1.5, color='#0066FF') +sns.lineplot(data=df_conc, x="timestamp", y="GLN_a", + errorbar=('sd'), label='[Gln]$_{a}$', lw=1.5, color="#44AB42F7") +sns.lineplot(data=df_conc, x="timestamp", y="GLU_a", + errorbar=('sd'), label='[Glu]$_{a}$', lw=1.5, color="dimgray") +sns.lineplot(data=df_conc, x="timestamp", y="GABA_i", + errorbar=('sd'), label='[GABA]$_{i}^c$', lw=1.5, color='orangered') + +# Format +plt.axhline(0, color='gray', lw=1, ls='-') +plt.xlabel('Time after bolus (min)') +plt.ylabel('Concentration (mM)', labelpad=3) +plt.xlim(0, 220) +plt.ylim([-0.5, 13]) +plt.xticks([0, 50, 100, 150, 200]) +leg = plt.legend(framealpha=1, fontsize=7, edgecolor="black", + bbox_to_anchor=(0.08, 0.9), loc="center left", + ncol=2, columnspacing=0.3, handlelength=1) +# plt.title("Pool concentrations (abs.) ", fontsize=10) +plt.grid(alpha=0.4, zorder=1) +leg.set_in_layout(False) + +# Validation with relative total concentrations (C) +# -------------- + +# Subplot +plt.subplot(1, 4, 3) + +# Glutamate +x = df_conc_norm["timestamp"] +y = df_conc_norm["GLU_tot"] +sns.lineplot(x=x, y=y, + errorbar=('sd'), label='[Glu]$_{tot}$ (sim.)', lw=1.5, color='#0066FF') + +# GABA +x = df_conc_norm["timestamp"] +y = df_conc_norm["GABA_tot"] +sns.lineplot(x=x, y=y, + errorbar=('sd'), label="[GABA]$_{tot}$ (sim.)", lw=1.5, color='orangered') + + +# Add experimental data +plt.errorbar(t_meas, exp_means['GLU'], xerr=t_meas_error, yerr=exp_errors["GLU"], + label="[Glu]$_{tot}$ (MRS)", fmt='o', lw=1.5, color="#001AFF", + capsize=2.5, capthick=1.5, markersize=2) +plt.errorbar(t_meas, exp_means['GABA'], xerr=t_meas_error, yerr=exp_errors["GABA"], + label="[GABA]$_{tot}$ (MRS)", fmt='o', lw=1.5, color="#BF0000", + capsize=2.5, capthick=1.5, markersize=2) +plt.axvline(t_meas, color='gray', lw=1, ls='--', zorder=-1) + +# Formatting +plt.axhline(0, color='gray', lw=1, ls='-') +plt.xlabel('Time after bolus (min)') +plt.ylabel(r"$\Delta$" + 'concentration (%)', labelpad=-5) +plt.xticks([0, 50, 100, 150, 200]) +plt.xlim(0, 220) +plt.ylim([-110, 13]) +leg = plt.legend(framealpha=1, edgecolor="black", fontsize=7, + bbox_to_anchor=(0.01, 0.16), loc="center left") +plt.grid(alpha=0.4) +leg.set_in_layout(False) + +# Relative fluxes (D) +# -------------- + +# Subplot +plt.subplot(1, 4, 4) + +# Plot +sns.lineplot(data=df_flux_norm, x="timestamp", y="GLU synaptic flux", errorbar=('sd'), + label='V$_{syn}^e$', lw=1.5, color="#0066FF") +sns.lineplot(data=df_flux_norm, x="timestamp", y="GABA synaptic flux", errorbar=('sd'), + label='V$_{syn}^i$', lw=1.5, color="#D60000F7") + + +# Format +plt.axhline(0, color='gray', lw=1, ls='-') +plt.xlabel('Time after bolus (min)') +plt.ylabel(r"$\Delta$" +'flux (%)', labelpad=-5) +plt.xticks([0, 50, 100, 150, 200]) +plt.xlim(0, 220) +plt.ylim([-15, 1.5]) +plt.grid(alpha=0.4) +leg = plt.legend(framealpha=1, fontsize=7, edgecolor="black", + bbox_to_anchor=(0.15, 0.07), loc="center left", + ncol=2, columnspacing=0.3, handlelength=1) +leg.set_in_layout(False) + +# Save +plt.tight_layout(w_pad=0.2) +plt.savefig(OUTDIR + f"figs/fig_timecourses_conf.pdf", + dpi=300, transparent=True) + +# %% +# =================================================================== +# SI figure: time-series +# =================================================================== + +# Figure +plt.figure(figsize=(8.25, 3.3)) + +# Fluxes +# ----------- + +# Subplot +plt.subplot(1, 3, 1) + +# Plot all fluxes +sns.lineplot(data=df_flux, x="timestamp", y="CMRglc", errorbar=('sd'), + label='CMR$_{glc}$', lw=1.5, color="#000000") +sns.lineplot(data=df_flux, x="timestamp", y="Ketone utilization", errorbar=('sd'), + label='CMR$_{Ket}$', lw=1.5, color="mediumvioletred") +sns.lineplot(data=df_flux, x="timestamp", y="PMAS in e", errorbar=('sd'), + label='PMAS$_e$', lw=1.5, color='navy') +sns.lineplot(data=df_flux, x="timestamp", y="PMAS in i", errorbar=('sd'), + label='PMAS$_i$', lw=1.5, color="maroon") +sns.lineplot(data=df_flux, x="timestamp", y="GLU synaptic flux", errorbar=('sd'), + label='V$_{syn}^e$', lw=1.5, color="#0066FF") +sns.lineplot(data=df_flux, x="timestamp", y="GABA synaptic flux", errorbar=('sd'), + label='V$_{syn}^i$', lw=1.5, color="#D60000F7") +sns.lineplot(data=df_flux, x="timestamp", y="GDH", errorbar=('sd'), + label='V$_{GDH}$', lw=1.5, color="#006B26FF") +sns.lineplot(data=df_flux, x="timestamp", y="GS", errorbar=('sd'), + label='V$_{GS}$', lw=1.5, color='#44AB42F7') +sns.lineplot(data=df_flux, x="timestamp", y="GAD65", errorbar=('sd'), + label='V$_{GAD65}$', lw=1.5, color="#DA6D85F7") +sns.lineplot(data=df_flux, x="timestamp", y="GAD67", errorbar=('sd'), + label='V$_{GAD67}$', lw=1.5, color="#D66DDAF7") +sns.lineplot(data=df_flux, x="timestamp", y="GT", errorbar=('sd'), + label='V$_{GT}$', lw=1.5, color="#4286ABF7") + + +# Format +plt.axhline(0, color='gray', lw=1, ls='-') +plt.xlabel('Time after bolus (min)') +plt.ylabel('Flux (mM/min)') +plt.xlim(0, 400) +plt.xticks([0, 100, 200, 300, 400]) +plt.grid(alpha=0.4) +leg = plt.legend(framealpha=1, fontsize=7, edgecolor="black", + bbox_to_anchor=(1.02, 0.42), loc="center left") +leg.set_in_layout(False) + +# Concentrations +# ----------- + +# Subplot +plt.subplot(1, 3, 2) + +# Plot all concentrations +sns.lineplot(data=df_conc, x="timestamp", y="GLU_e", + errorbar=('sd'), label='[Glu]$_{e}$', lw=1.5, color='#0066FF') +sns.lineplot(data=df_conc, x="timestamp", y="GLN_a", + errorbar=('sd'), label='[Gln]$_{a}$', lw=1.5, color="#44AB42F7") +sns.lineplot(data=df_conc, x="timestamp", y="GLU_a", + errorbar=('sd'), label='[Glu]$_{a}$', lw=1.5, color='dimgray') +sns.lineplot(data=df_conc, x="timestamp", y="GABA_i", + errorbar=('sd'), label='[GABA]$_{i}^c$', lw=1.5, color='orangered') +sns.lineplot(data=df_conc, x="timestamp", y="GABA_v", + errorbar=('sd'), label='[GABA]$_{i}^v$', lw=1.5, color='maroon') +sns.lineplot(data=df_conc, x="timestamp", y="GLU_i", + errorbar=('sd'), label='[Glu]$_{i}$', lw=1.5, color='dodgerblue') +sns.lineplot(data=df_conc, x="timestamp", y="GLN_e", + errorbar=('sd'), label='[Gln]$_{e}$', lw=1.5, color='yellowgreen') +sns.lineplot(data=df_conc, x="timestamp", y="GLN_i", + errorbar=('sd'), label='[Gln]$_{i}$', lw=1.5, color='lightcoral') +sns.lineplot(data=df_conc, x="timestamp", y="KET_b", + errorbar=('sd'), label='[KetB]$_{e}$', lw=1.5, color='darkorange') + + +# Format +plt.axhline(0, color='gray', lw=1, ls='-') +plt.xlabel('Time after bolus (min)') +plt.ylabel('Concentration (mM)') +plt.xlim(0, 400) +plt.xticks([0, 100, 200, 300, 400]) +leg = plt.legend(framealpha=1, fontsize=7, edgecolor="black", + bbox_to_anchor=(1.02, 0.42), loc="center left") +plt.grid(alpha=0.4, zorder=1) +leg.set_in_layout(False) + +# Subplot +plt.subplot(1, 3, 3) + +# Plot all concentrations +sns.lineplot(data=df_conc, x="timestamp", y="GABA_i", + errorbar=('sd'), label='[GABA]$_{i}^c$', lw=1.5, color='orangered') +sns.lineplot(data=df_conc, x="timestamp", y="GABA_v", + errorbar=('sd'), label='[GABA]$_{i}^v$', lw=1.5, color='maroon') +sns.lineplot(data=df_conc, x="timestamp", y="GLU_a", + errorbar=('sd'), label='[Glu]$_{a}$', lw=1.5, color='dimgray') +sns.lineplot(data=df_conc, x="timestamp", y="GLU_i", + errorbar=('sd'), label='[Glu]$_{i}$', lw=1.5, color='dodgerblue') +sns.lineplot(data=df_conc, x="timestamp", y="GLN_e", + errorbar=('sd'), label='[Gln]$_{e}$', lw=1.5, color='yellowgreen') +sns.lineplot(data=df_conc, x="timestamp", y="GLN_i", + errorbar=('sd'), label='[Gln]$_{i}$', lw=1.5, color='lightcoral') +sns.lineplot(data=df_conc, x="timestamp", y="KET_b", + errorbar=('sd'), label='[KetB]$_{e}$', lw=1.5, color='darkorange') + +# Format +plt.axhline(0, color='gray', lw=1, ls='-') +plt.xlabel('Time after bolus (min)') +plt.ylabel('Concentration (mM)') +plt.xlim(0, 400) +plt.xticks([0, 100, 200, 300, 400]) +leg = plt.legend(framealpha=1, fontsize=7, edgecolor="black", + bbox_to_anchor=(0.42, 0.72), loc="center left") +plt.grid(alpha=0.4, zorder=1) +leg.set_in_layout(False) + +# Save +plt.tight_layout(w_pad=5.5) +plt.savefig(OUTDIR + f"figs/SI_fig_timecourses_all.pdf", dpi=300, transparent=True) + + +# %% +# =================================================================== +# SI figure: times-series (scaled to compartmental concentrations) +# =================================================================== + +# Scaling factors +vol_e = 0.56 +vol_i = 0.14 +vol_a = 0.2 + +# Figure +plt.figure(figsize=(8.25, 3.3)) + +# Fluxes +# ----------- + +# Subplot +plt.subplot(1, 3, 1) + +# Plot all fluxes +sns.lineplot(x=df_flux["timestamp"], y=df_flux["CMRglc"], errorbar=('sd'), + label='CMR$_{glc}$', lw=1.5, color="#000000") +sns.lineplot(x=df_flux["timestamp"], y=df_flux["Ketone utilization"], errorbar=('sd'), + label='CMR$_{Ket}$', lw=1.5, color="mediumvioletred") +sns.lineplot(x=df_flux["timestamp"], y=df_flux["PMAS in e"]/vol_e, errorbar=('sd'), + label='PMAS$_e$', lw=1.5, color='navy') +sns.lineplot(x=df_flux["timestamp"], y=df_flux["PMAS in i"]/vol_i, errorbar=('sd'), + label='PMAS$_i$', lw=1.5, color="maroon") +sns.lineplot(x=df_flux["timestamp"], y=df_flux["GLU synaptic flux"]/vol_e, errorbar=('sd'), + label='V$_{syn}^e$', lw=1.5, color="#0066FF") +sns.lineplot(x=df_flux["timestamp"], y=df_flux["GABA synaptic flux"]/vol_i, errorbar=('sd'), + label='V$_{syn}^i$', lw=1.5, color="#D60000F7") +sns.lineplot(x=df_flux["timestamp"], y=df_flux["GDH"]/vol_a, errorbar=('sd'), + label='V$_{GDH}$', lw=1.5, color="#006B26FF") +sns.lineplot(x=df_flux["timestamp"], y=df_flux["GS"]/vol_a, errorbar=('sd'), + label='V$_{GS}$', lw=1.5, color='#44AB42F7') +sns.lineplot(x=df_flux["timestamp"], y=df_flux["GAD65"]/vol_i, errorbar=('sd'), + label='V$_{GAD65}$', lw=1.5, color="#DA6D85F7") +sns.lineplot(x=df_flux["timestamp"], y=df_flux["GAD67"]/vol_i, errorbar=('sd'), + label='V$_{GAD67}$', lw=1.5, color="#D66DDAF7") +sns.lineplot(x=df_flux["timestamp"], y=df_flux["GT"]/vol_i, errorbar=('sd'), + label='V$_{GT}$', lw=1.5, color="#4286ABF7") + + +# Format +plt.axhline(0, color='gray', lw=1, ls='-') +plt.xlabel('Time after bolus (min)') +plt.ylabel('Flux (mM/min)') +plt.xlim(0, 400) +plt.xticks([0, 100, 200, 300, 400]) +plt.grid(alpha=0.4) +leg = plt.legend(framealpha=1, fontsize=7, edgecolor="black", + bbox_to_anchor=(1.02, 0.42), loc="center left") +leg.set_in_layout(False) + +# Concentrations +# ----------- + +# Subplot +plt.subplot(1, 3, 2) + +# Plot all concentrations +sns.lineplot(x=df_conc["timestamp"], y=df_conc["GLU_e"]/vol_e, + errorbar=('sd'), label='[Glu]$_{e}$', lw=1.5, color='#0066FF') +sns.lineplot(x=df_conc["timestamp"], y=df_conc["GLN_a"]/vol_a, + errorbar=('sd'), label='[Gln]$_{a}$', lw=1.5, color="#44AB42F7") +sns.lineplot(x=df_conc["timestamp"], y=df_conc["GLU_a"]/vol_a, + errorbar=('sd'), label='[Glu]$_{a}$', lw=1.5, color='dimgray') +sns.lineplot(x=df_conc["timestamp"], y=df_conc["GABA_i"]/vol_i, + errorbar=('sd'), label='[GABA]$_{i}^c$', lw=1.5, color='orangered') +sns.lineplot(data=df_conc, x="timestamp", y="GABA_v", + errorbar=('sd'), label='[GABA]$_{i}^v$', lw=1.5, color='maroon') +sns.lineplot(x=df_conc["timestamp"], y=df_conc["GLU_i"]/vol_i, + errorbar=('sd'), label='[Glu]$_{i}$', lw=1.5, color='dodgerblue') +sns.lineplot(x=df_conc["timestamp"], y=df_conc["GLN_e"]/vol_e, + errorbar=('sd'), label='[Gln]$_{e}$', lw=1.5, color='yellowgreen') +sns.lineplot(x=df_conc["timestamp"], y=df_conc["GLN_i"]/vol_i, + errorbar=('sd'), label='[Gln]$_{i}$', lw=1.5, color='lightcoral') +sns.lineplot(x=df_conc["timestamp"], y=df_conc["KET_b"], + errorbar=('sd'), label='[KetB]$_{e}$', lw=1.5, color='darkorange') + +# Format +plt.axhline(0, color='gray', lw=1, ls='-') +plt.xlabel('Time after bolus (min)') +plt.ylabel('Concentration (mM)') +plt.xlim(0, 400) +plt.xticks([0, 100, 200, 300, 400]) +leg = plt.legend(framealpha=1, fontsize=7, edgecolor="black", + bbox_to_anchor=(1.02, 0.42), loc="center left") +plt.grid(alpha=0.4, zorder=1) +leg.set_in_layout(False) + +# Subplot +plt.subplot(1, 3, 3) + +# Plot all concentrations +sns.lineplot(x=df_conc["timestamp"], y=df_conc["GABA_i"]/vol_i, + errorbar=('sd'), label='[GABA]$_{i}^c$', lw=1.5, color='orangered') +sns.lineplot(data=df_conc, x="timestamp", y="GABA_v", + errorbar=('sd'), label='[GABA]$_{i}^v$', lw=1.5, color='maroon') +sns.lineplot(x=df_conc["timestamp"], y=df_conc["GLU_a"]/vol_a, + errorbar=('sd'), label='[Glu]$_{a}$', lw=1.5, color='dimgray') +sns.lineplot(x=df_conc["timestamp"], y=df_conc["GLU_i"]/vol_i, + errorbar=('sd'), label='[Glu]$_{i}$', lw=1.5, color='dodgerblue') +sns.lineplot(x=df_conc["timestamp"], y=df_conc["GLN_e"]/vol_e, + errorbar=('sd'), label='[Gln]$_{e}$', lw=1.5, color='yellowgreen') +sns.lineplot(x=df_conc["timestamp"], y=df_conc["GLN_i"]/vol_i, + errorbar=('sd'), label='[Gln]$_{i}$', lw=1.5, color='lightcoral') +sns.lineplot(x=df_conc["timestamp"], y=df_conc["KET_b"], + errorbar=('sd'), label='[KetB]$_{e}$', lw=1.5, color='darkorange') + +# Format +plt.axhline(0, color='gray', lw=1, ls='-') +plt.xlabel('Time after bolus (min)') +plt.ylabel('Concentration (mM)') +plt.xlim(0, 400) +plt.xticks([0, 100, 200, 300, 400]) +leg = plt.legend(framealpha=1, fontsize=7, edgecolor="black", + bbox_to_anchor=(0.42, 0.72), loc="center left") +plt.grid(alpha=0.4, zorder=1) +leg.set_in_layout(False) + +# Save +plt.tight_layout(w_pad=5.5) +plt.savefig(OUTDIR + f"figs/SI_fig_timecourses_all_compartmental.pdf", + dpi=300, transparent=True) + +# %% +# =================================================================== +# Figure 4: MCA +# =================================================================== + +# Open MCA results +data_mca_conc_ctrl = pd.read_csv(OUTDIR + f"model_mca_concentration.csv"); + +# Settings +vmin = -10.1 +vmax = 10.1 + +# Figure +plt.figure(figsize=(8.25, 4.)) + +# Transform data for plotting +df_plot = data_mca_conc_ctrl \ + .set_index("state") \ + .loc[ + ["GLU_e(t)", "GLU_i(t)", "GLU_a(t)", "GLN_e(t)", + "GLN_i(t)","GLN_a(t)", "GABA_i(t)", "GABA_v(t)"], + ["c_vpmas", + "Vmax_GAD65", "Vmax_GAD67", "Vmax_GT", "Vmax_GDH", "Vmax_GS", "KM_GAD65", + "KM_GAD67", "KM_GT", "KM_GDH", "h_GDH", "KM_GS", "KI_GS", "k_syn_e", "k_syn_i"] + ] + +# Annot data +annot_ctrl = np.where(~np.isclose(df_plot, 0, atol=2e-1), df_plot.round(1).astype(str), "") + +# Plot heatmap +sns.heatmap(df_plot, annot=annot_ctrl, fmt="", cmap="seismic", square=True, + vmin=vmin, vmax=vmax, annot_kws={"size": 7}, + cbar_kws={"label": "Concentration control coefficient", "shrink": 1, + "ticks": np.arange(vmin+0.1, vmax+0.9, 5)}) + +# Format +xticks = ["V$_{PMAS}$", "V$_{max}^{GAD65}$", "V$_{max}^{GAD67}$", + "V$_{max}^{GT}$", "V$_{max}^{GDH}$", "V$_{max}^{GS}$", "K$_M^{GAD65}$", "K$_M^{GAD67}$", "K$_M^{GT}$", + "K$_M^{GDH}$", "h$_{GDH}$", "K$_M^{GS}$", "K$_I^{GS}$", + "k$_{syn}^{e}$", "k$_{syn}^{i}$"] +yticks = ["\n\nGLU$_{e}$", "\n\nGLU$_{i}$", "\n\nGLU$_{a}$", "\n\nGLN$_{e}$", + "\n\nGLN$_{i}$","\n\nGLN$_{a}$", "\n\nGABA$_{i}^{c}$", "\n\nGABA$_{i}^{v}$"] + +# plt.title("Concentration control") +plt.xticks(np.arange(len(xticks) + 1), xticks + [None], rotation=60, ha="left") +plt.yticks(np.arange(len(yticks) + 1), yticks + [None], rotation=0, va="center") + +plt.xlabel("Parameter") +plt.ylabel("Metabolite pool") +plt.grid(alpha=0.5) + +for spine in plt.gca().spines.values(): + spine.set_visible(True) + spine.set_linewidth(0.5) + spine.set_color("black") + +for spine in plt.gca().spines.values(): + spine.set_visible(True) + spine.set_linewidth(0.5) + spine.set_color("black") + +# Save +plt.tight_layout(rect=[0.01, -0.0, 1.02, 1.0], w_pad=-20) +plt.savefig(OUTDIR + f"figs/fig_mca_concentration.pdf", dpi=300, transparent=True) diff --git a/reproduce_paper/code/2_plot_results.py b/reproduce_paper/code/2_plot_results.py deleted file mode 100644 index 8bce708..0000000 --- a/reproduce_paper/code/2_plot_results.py +++ /dev/null @@ -1,470 +0,0 @@ -""" - -@author: Botond B. Antal - -January 2026 - -This script plots the simulated results. - -""" - -import os -import pandas as pd -import numpy as np -import matplotlib.pyplot as plt -import seaborn as sns - -# Directories -OUTDIR = "" + "results/" #TODO: insert your filepath here before results/ - -# rcParams -plt.rcParams.update({'font.size': 10}) - -# %% -# =================================================================== -# Figure 1: ketones -# =================================================================== - -# Load data -data_bloodbhb_datapoints = pd.read_csv(OUTDIR + "martin-arrowsmith_etal_blood_bhb_timecourse.csv") -data_bloodbhb_timeseries = pd.read_csv(OUTDIR + f"model_blood_bhb_input_timeseries.csv") -data_brainbhb_datapoints = pd.read_csv(OUTDIR + "mujica-parodi_etal_brain_bhb_timeseries.csv") -data_brainbhb_timeseries = pd.read_csv(OUTDIR + f"model_timeseries.csv")[["timestamp", "KET_b"]].set_axis(["time", "concentration"], axis=1) - -# Figure -plt.figure(figsize=(8.25, 3.5)) - -# Subplot 1 -plt.subplot(1, 3, 1) -plt.scatter( - data_bloodbhb_datapoints['time'], - data_bloodbhb_datapoints['concentration'], - label='Measured\n(Martin-\nArrowsmith\net al.)', - color='black', s=80, zorder=2) -plt.plot( - data_bloodbhb_timeseries['time'], - data_bloodbhb_timeseries['concentration'], - linewidth=2, label='Fitted (double\nexponential)', - color='crimson', zorder=3) -plt.xlabel('Time after bolus (min)') -plt.ylabel('Concentration (mM)') -plt.title(r"Blood D-$\beta$HB") -plt.xlim(0, 220) -plt.xticks([0, 50, 100, 150, 200]) -leg = plt.legend(framealpha=1, fontsize=8, bbox_to_anchor=(0.08, 0.2), - edgecolor="black", loc="center left") -plt.grid(alpha=0.4, zorder=1) -leg.set_in_layout(False) - -# Subplot 2 -plt.subplot(1, 3, 2) -plt.scatter( - data_brainbhb_datapoints['x'], - data_brainbhb_datapoints['y']/12, - label='Measured and rescaled\n(Mujica-Parodi et al.)', - color='black', s=80, zorder=2) -plt.plot( - data_brainbhb_timeseries['time'], - data_brainbhb_timeseries['concentration'], - linewidth=2, label='Simulated', color='dodgerblue', zorder=3) -plt.xlabel('Time after bolus (min)') -plt.ylabel('Concentration (mM)') -plt.title(r"Brain D-$\beta$HB") -plt.xlim(0, 220) -plt.xticks([0, 50, 100, 150, 200]) -plt.legend(framealpha=1, fontsize=8, edgecolor="black") -plt.grid(alpha=0.4, zorder=1) - -# Subplot 3 -plt.subplot(1, 3, 3) -plt.plot( - data_bloodbhb_timeseries['time'], - data_bloodbhb_timeseries['concentration']/data_bloodbhb_timeseries['concentration'].max(), - linewidth=2, label='Blood', color='crimson', zorder=2) -plt.plot( - data_brainbhb_timeseries['time'], - data_brainbhb_timeseries['concentration']/data_brainbhb_timeseries['concentration'].max(), - linewidth=2, label='Brain', color='dodgerblue', zorder=3) -plt.xlabel('Time after bolus (min)') -plt.ylabel('Concentration (relative)') -plt.title("Normalized\n" + r'Blood & Brain D-$\beta$HB') -plt.xlim(0, 100) -plt.legend(framealpha=1, fontsize=8, edgecolor="black") -plt.grid(alpha=0.4, zorder=0) - -# Save -plt.tight_layout() -plt.savefig(OUTDIR + f"fig_ketone_pk.pdf", - dpi=300, transparent=True) - -# %% -# =================================================================== -# Figure 2: time-series and validation -# =================================================================== - -# Load simulation results -data_concs = pd.read_csv(OUTDIR + f"model_timeseries.csv") -data_fluxes = pd.read_csv(OUTDIR + f"model_fluxes.csv") - -# Experiment data -exp_means = {"GLU": -10.6, "GABA": -36.9} -exp_errors = {"GLU": 1.2*1.96, "GABA": 4.2*1.96} # 95% CI! - -# Figure -plt.figure(figsize=(8.25, 3.3)) - -# Fluxes -# -------------- - -# Subplot -plt.subplot(1, 3, 1) - -# Plot -plt.plot(data_fluxes["timestamp"], data_fluxes["CMRglc"], - label='CMR$_{glc}$', lw=2, color="#000000") -plt.plot(data_fluxes["timestamp"], data_fluxes["GS"], - label='V$_{GS}$', lw=2, linestyle='-', color="#44AB42F7") -plt.plot(data_fluxes["timestamp"], data_fluxes["PMAS in e"], - label='V$_{PMAS}^e$', lw=2, color='#0066FF') -plt.plot(data_fluxes["timestamp"], data_fluxes["PMAS in i"], - label='V$_{PMAS}^i$', lw=2, color="#D60000F7") -plt.plot(data_fluxes["timestamp"], data_fluxes["Ketone utilization"], - label='CMR$_{Ket}$', lw=2, linestyle='--', color="#686868") -# plt.plot(data_fluxes["timestamp"], data_fluxes["GAD65"], -# label='V$_{GAD65}$', lw=2, color="orange") -# plt.plot(data_fluxes["timestamp"], data_fluxes["GAD67"], -# label='V$_{GAD67}$', lw=2, color="indianred") -# plt.plot(data_fluxes["timestamp"], data_fluxes["GT"], -# label='V$_{GT}$', lw=2, linestyle='--', color="#B51861F7") - -# Format -plt.xlabel('Time after bolus (min)') -plt.ylabel('Flux (mM/min)') -plt.xticks([0, 50, 100, 150, 200]) -plt.xlim(0, 220) -plt.grid(alpha=0.4) -plt.title("Fluxes") -leg = plt.legend(framealpha=1, fontsize=7, edgecolor="black", - bbox_to_anchor=(0.52, 0.48), loc="center left") -leg.set_in_layout - -# Print values -print(f"CMRglc starts from {data_fluxes['CMRglc'][0]:.3f} mM/min " - f"and is lowest at {data_fluxes['CMRglc'].min():.3f} mM/min " - f"which occurs at t={data_fluxes['timestamp'][data_fluxes['CMRglc'].idxmin()]:.2f}) min") -print(f"CMRket starts from {data_fluxes['Ketone utilization'][0]:.3f} mM/min " - f"and is highest at {data_fluxes['Ketone utilization'].max():.3f} mM/min " - f"which occurs at t={data_fluxes['timestamp'][data_fluxes['Ketone utilization'].idxmax()]:.2f} min") -print(f"PMAS_e starts from {data_fluxes['PMAS in e'][0]:.3f} mM/min " - f"and is lowest at {data_fluxes['PMAS in e'].min():.3f} mM/min " - f"which occurs at t={data_fluxes['timestamp'][data_fluxes['PMAS in e'].idxmin()]:.2f} min") -print(f"PMAS_i starts from {data_fluxes['PMAS in i'][0]:.3f} mM/min " - f"and is lowest at {data_fluxes['PMAS in i'].min():.3f} mM/min " - f"which occurs at t={data_fluxes['timestamp'][data_fluxes['PMAS in i'].idxmin()]:.2f} min") -print(f"VGDH starts from {data_fluxes['GDH'][0]:.3f} mM/min " - f"and is lowest at {data_fluxes['GDH'].min():.3f} mM/min " - f"which occurs at t={data_fluxes['timestamp'][data_fluxes['GDH'].idxmin()]:.2f} min") -print(f"VGS starts from {data_fluxes['GS'][0]:.3f} mM/min " - f"and is lowest at {data_fluxes['GS'].min():.3f} mM/min " - f"which occurs at t={data_fluxes['timestamp'][data_fluxes['GS'].idxmin()]:.2f} min") -print(f"VGAD65 starts from {data_fluxes['GAD65'][0]:.3f} mM/min " - f"and is lowest at {data_fluxes['GAD65'].min():.3f} mM/min " - f"which occurs at t={data_fluxes['timestamp'][data_fluxes['GAD65'].idxmin()]:.2f} min") -print(f"VGAD67 starts from {data_fluxes['GAD67'][0]:.3f} mM/min " - f"and is lowest at {data_fluxes['GAD67'].min():.3f} mM/min " - f"which occurs at t={data_fluxes['timestamp'][data_fluxes['GAD67'].idxmin()]:.2f} min") -print(f"VGT starts from {data_fluxes['GT'][0]:.3f} mM/min " - f"and is lowest at {data_fluxes['GT'].min():.3f} mM/min " - f"which occurs at t={data_fluxes['timestamp'][data_fluxes['GT'].idxmin()]:.2f} min") - -# Compartmental concentrations -# -------------- - -# Subplot -plt.subplot(1, 3, 2) - -# Plot -plt.plot(data_concs["timestamp"], data_concs["GLU_e"], - label='[Glu]$_{e}$', lw=2, color='#0066FF') -plt.plot(data_concs["timestamp"], data_concs["GLN_a"], - label='[Gln]$_{a}$', lw=2, color="#44AB42F7") -plt.plot(data_concs["timestamp"], data_concs["GLU_a"], - label='[Glu]$_{a}$', lw=2, color="dimgray") -plt.plot(data_concs["timestamp"], data_concs["GABA_i"], - label='[GABA]$_{i}^c$', lw=2, color='orangered') -plt.plot(data_concs["timestamp"], data_concs["GABA_v"], - label='[GABA]$_{i}^v$', lw=2, color='maroon') - -# Format -plt.xlabel('Time after bolus (min)') -plt.ylabel('Concentration (mM)') -plt.xlim(0, 220) -plt.xticks([0, 50, 100, 150, 200]) -leg = plt.legend(framealpha=1, fontsize=7, edgecolor="black", - bbox_to_anchor=(0.56, 0.62), loc="center left") -plt.title("Pool concentrations (abs.) ") -plt.grid(alpha=0.4, zorder=1) -leg.set_in_layout(False) - -# Print values -print(f"GLU_e starts from {data_concs['GLU_e'][0]:.3f} mM " - f"and is lowest at {data_concs['GLU_e'].min():.3f} mM " - f"(-{100 * (data_concs['GLU_e'][0] - data_concs['GLU_e'].min()) / data_concs['GLU_e'][0]:.1f}%) " - f"which occurs at t={data_concs['timestamp'][data_concs['GLU_e'].idxmin()]:.2f} min") -print(f"GLU_a starts from {data_concs['GLU_a'][0]:.3f} mM " - f"and is lowest at {data_concs['GLU_a'].min():.3f} mM " - f"(-{100 * (data_concs['GLU_a'][0] - data_concs['GLU_a'].min()) / data_concs['GLU_a'][0]:.1f}%) " - f"which occurs at t={data_concs['timestamp'][data_concs['GLU_a'].idxmin()]:.2f} min") -print(f"GLU_a starts from {data_concs['GLU_a'][0]:.3f} mM " - f"and is highest at {data_concs['GLU_a'].max():.3f} mM " - f"(-{100 * (data_concs['GLU_a'][0] - data_concs['GLU_a'].max()) / data_concs['GLU_a'][0]:.1f}%) " - f"which occurs at t={data_concs['timestamp'][data_concs['GLU_a'].idxmax()]:.2f} min") -print(f"GLN_a starts from {data_concs['GLN_a'][0]:.3f} mM " - f"and is highest at {data_concs['GLN_a'].max():.3f} mM " - f"(-{100 * (data_concs['GLN_a'][0] - data_concs['GLN_a'].max()) / data_concs['GLN_a'][0]:.1f}%) " - f"which occurs at t={data_concs['timestamp'][data_concs['GLN_a'].idxmax()]:.2f} min") -print(f"GABA_i starts from {data_concs['GABA_i'][0]:.3f} mM " - f"and is lowest at {data_concs['GABA_i'].min():.3f} mM " - f"(-{100 * (data_concs['GABA_i'][0] - data_concs['GABA_i'].min()) / data_concs['GABA_i'][0]:.1f}%) " - f"which occurs at t={data_concs['timestamp'][data_concs['GABA_i'].idxmin()]:.2f} min") -print(f"GABA_v starts from {data_concs['GABA_v'][0]:.3f} mM " - f"and is lowest at {data_concs['GABA_v'].min():.3f} mM " - f"(-{100 * (data_concs['GABA_v'][0] - data_concs['GABA_v'].min()) / data_concs['GABA_v'][0]:.1f}%) " - f"which occurs at t={data_concs['timestamp'][data_concs['GABA_v'].idxmin()]:.2f} min") - - -# Validation with total concentrations -# -------------- - -# Subplot -plt.subplot(1, 3, 3) - -# Glutamate -cols = ["GLU_e", "GLU_i", "GLU_a"] -x = data_concs["timestamp"] -y = 100 * data_concs[cols].sum(axis=1) / \ - data_concs.loc[0, cols].sum() - 100 -plt.plot(x, y, label='[Glu]$_{tot}$ (sim.)', lw=2, color='#0066FF') -y_60_glu = y[(data_concs['timestamp'] - 60).abs().idxmin()] - -# GABA -cols = ["GABA_i", "GABA_v"] -x = data_concs["timestamp"] -y = 100 * data_concs[cols].sum(axis=1) / \ - data_concs.loc[0, cols].sum() - 100 -plt.plot(x, y, label="[GABA]$_{tot}$ (sim.)", lw=2, color='orangered') -y_60_gaba = y[(data_concs['timestamp'] - 60).abs().idxmin()] - -# Add experimental data -plt.errorbar(60, exp_means['GLU'], yerr=exp_errors["GLU"], - label="[Glu]$_{tot}$ (meas.)", fmt='o', lw=2, color="#001AFF", - capsize=5, capthick=2) -plt.errorbar(60, exp_means['GABA'], yerr=exp_errors["GABA"], - label="[GABA]$_{tot}$ (meas.)", fmt='o', lw=2, color="#BF0000", - capsize=5, capthick=2) -plt.axvline(60, color='black', lw=1.5, ls='--', zorder=-1) - -# Formatting -plt.axhline(0, color='black', lw=1, ls='--') -plt.xlabel('Time after bolus (min)') -plt.ylabel(r"$\Delta$%" + ' concentration (%)') -plt.xticks([0, 50, 100, 150, 200]) -plt.xlim(0, 220) -plt.ylim(-100, 13) -leg = plt.legend(framealpha=1, edgecolor="black", fontsize=7, - bbox_to_anchor=(0.01, 0.19), loc="center left") -plt.grid(alpha=0.4) -plt.title("Total concentrations (" + r'$\Delta$%)') -leg.set_in_layout(False) - -# Save -# plt.subplots_adjust(wspace=0.9) -plt.tight_layout(w_pad=1) -plt.savefig(OUTDIR + f"fig_timecourses.pdf", - dpi=300, transparent=True) - -# Print values -print(f"Glutamate: {y_60_glu:.1f}%") -print(f"GABA: {y_60_gaba:.1f}%") - -# %% -# =================================================================== -# SI figure: time-series -# =================================================================== - -# Load simulation results -data_concs = pd.read_csv(OUTDIR + f"model_timeseries.csv") -data_fluxes = pd.read_csv(OUTDIR + f"model_fluxes.csv") - -# Figure -plt.figure(figsize=(8.25, 3.3)) - -# Fluxes -# ----------- - -# Subplot -plt.subplot(1, 3, 1) - -# Plot all fluxes -plt.plot(data_fluxes["timestamp"], (data_fluxes["CMRglc"]), - label='CMR$_{glc}$', lw=2, color="#424447") -plt.plot(data_fluxes["timestamp"], data_fluxes["Ketone utilization"], - label='CMR$_{Ket}$', lw=2, color="#A22727") -plt.plot(data_fluxes["timestamp"], data_fluxes["PMAS in e"], - label='PMAS$_e$', lw=2, color='#001DD6F7') -plt.plot(data_fluxes["timestamp"], data_fluxes["PMAS in i"], - label='PMAS$_i$', lw=2, color="#D60000F7") -plt.plot(data_fluxes["timestamp"], data_fluxes["GLU cycling"], - label='V$_{syn}^e$', lw=2, color="#00BFFF") -plt.plot(data_fluxes["timestamp"], data_fluxes["GABA cycling"], - label='V$_{syn}^i$', lw=2, color="#FF00F7") -plt.plot(data_fluxes["timestamp"], data_fluxes["GDH"], - label='V$_{GDH}$', lw=2, color="#006B26FF") -plt.plot(data_fluxes["timestamp"], data_fluxes["GS"], - label='V$_{GS}$', lw=2, color='#44AB42F7') -plt.plot(data_fluxes["timestamp"], data_fluxes["GAD65"], - label='V$_{GAD65}$', lw=2, color="#DA6D85F7") -plt.plot(data_fluxes["timestamp"], data_fluxes["GAD67"], - label='V$_{GAD67}$', lw=2, color="#D66DDAF7") -plt.plot(data_fluxes["timestamp"], data_fluxes["GT"], - label='V$_{GT}$', lw=2, color="#4286ABF7") - - -# Format -plt.xlabel('Time after bolus (min)') -plt.ylabel('Flux (mM/min)') -plt.xlim(0, 400) -plt.grid(alpha=0.4) -plt.title("Fluxes") -leg = plt.legend(framealpha=1, fontsize=7, edgecolor="black", - bbox_to_anchor=(1.02, 0.42), loc="center left") -leg.set_in_layout(False) - -# Concentrations -# ----------- - -# Subplot -plt.subplot(1, 3, 2) - -# Plot all concentrations -plt.plot(data_concs["timestamp"], data_concs["GLU_e"], - label='[Glu]$_{e}$', lw=2, color='#0066FF') -plt.plot(data_concs["timestamp"], data_concs["GLN_a"], - label='[Gln]$_{a}$', lw=2, color="#44AB42F7") -plt.plot(data_concs["timestamp"], data_concs["GABA_i"], - label='[GABA]$_{i}$', lw=2, color='orangered') -plt.plot(data_concs["timestamp"], data_concs["GABA_v"], - label='[GABA]$_{v}$', lw=2, color='maroon') -plt.plot(data_concs["timestamp"], data_concs["GLU_a"], - label='[Glu]$_{a}$', lw=2, color='deepskyblue') -plt.plot(data_concs["timestamp"], data_concs["GLU_i"], - label='[Glu]$_{i}$', lw=2, color='dodgerblue') -plt.plot(data_concs["timestamp"], data_concs["GLN_e"], - label='[Gln]$_{e}$', lw=2, color='yellowgreen') -plt.plot(data_concs["timestamp"], data_concs["GLN_i"], - label='[Gln]$_{i}$', lw=2, color='lightcoral') -# plt.plot(data_concs["timestamp"], data_concs["GLU_im"], -# label='[Glu]$_{im}$', lw=2, color='blueviolet') -plt.plot(data_concs["timestamp"], data_concs["KET_b"], - label='[KetB]$_{e}$', lw=2, color='darkorange') - -# Format -plt.xlabel('Time after bolus (min)') -plt.ylabel('Concentration (mM)') -plt.xlim(0, 400) -# plt.xticks([0, 50, 100, 150, 200]) -leg = plt.legend(framealpha=1, fontsize=7, edgecolor="black", - bbox_to_anchor=(1.02, 0.42), loc="center left") -plt.title("Pool concentrations") -plt.grid(alpha=0.4, zorder=1) -leg.set_in_layout(False) - - -# Subplot -plt.subplot(1, 3, 3) - -# Plot all concentrations -plt.plot(data_concs["timestamp"], data_concs["GABA_i"], - label='[GABA]$_{i}$', lw=2, color='orangered') -plt.plot(data_concs["timestamp"], data_concs["GLU_a"], - label='[Glu]$_{a}$', lw=2, color='deepskyblue') -plt.plot(data_concs["timestamp"], data_concs["GLN_i"], - label='[Gln]$_{i}$', lw=2, color='lightcoral') -plt.plot(data_concs["timestamp"], data_concs["GLU_i"], - label='[Glu]$_{i}$', lw=2, color='dodgerblue') -plt.plot(data_concs["timestamp"], data_concs["GABA_v"], - label='[GABA]$_{v}$', lw=2, color='maroon') -plt.plot(data_concs["timestamp"], data_concs["KET_b"], - label='[KET]$_{b}$', lw=2, color='darkorange') - -# Format -plt.xlabel('Time after bolus (min)') -plt.ylabel('Concentration (mM)') -plt.xlim(0, 400) -# plt.xticks([0, 50, 100, 150, 200]) -leg = plt.legend(framealpha=1, fontsize=7, edgecolor="black", - bbox_to_anchor=(0.59, 0.72), loc="center left") -plt.title("(Small) pool concentrations ") -plt.grid(alpha=0.4, zorder=1) -leg.set_in_layout(False) - -# Save -plt.tight_layout(w_pad=5.5) -plt.savefig(OUTDIR + f"SI_fig_timecourses_all.pdf", dpi=300, transparent=True) - - -# %% -# =================================================================== -# Figure 3: MCA -# =================================================================== - -# Open MCA results -data_mca_conc_ctrl = pd.read_csv(OUTDIR + f"model_mca_concentration.csv"); - -# Figure -plt.figure(figsize=(8.25, 3.8)) - -# Concentration control -# ---------- - -# Transform -data_mca_conc_ctrl.set_index("state", inplace=True) # Set indexes -data_mca_conc_ctrl = data_mca_conc_ctrl \ - .loc[:, ["c_vpmas", "k_syn_e", "k_syn_i", "Vmax_GAD65", "Vmax_GAD67", "Vmax_GT", "Vmax_GDH", "Vmax_GS", "KM_GAD65", "KM_GAD67", "KM_GT", "KM_GDH", "KM_GS", "KI_GS"]] # Keep only relevant columns -data_mca_conc_ctrl = data_mca_conc_ctrl \ - .loc[["GLU_e(t)", "GLU_i(t)", "GLU_a(t)", "GLN_e(t)", "GLN_i(t)","GLN_a(t)", "GABA_i(t)", "GABA_v(t)"], :] # Keep only relevant rows - -# Subplot -# plt.subplot(1, 2, 1) - -# Annot data -annot_ctrl = np.where(~np.isclose(data_mca_conc_ctrl, 0, atol=2e-1), data_mca_conc_ctrl.round(1).astype(str), "") - -# Plot heatmap -sns.heatmap(data_mca_conc_ctrl, annot=annot_ctrl, fmt="", cmap="seismic", square=True, - vmin=-25, vmax=25, annot_kws={"size": 8}, - cbar_kws={"label": "Concentration control coefficient", "shrink": 1}) - -# Format -xticks = ["V$_{PMAS}$", "k$_{syn}^{e}$", "k$_{syn}^{i}$", "V$_{max}^{GAD65}$", "V$_{max}^{GAD67}$", - "V$_{max}^{GT}$", "V$_{max}^{GDH}$", "V$_{max}^{GS}$", "K$_M^{GAD65}$", "K$_M^{GAD67}$", "K$_M^{GT}$", - "K$_M^{GDH}$", "K$_M^{GS}$", "K$_I^{GS}$", ] -yticks = ["\n\nGLU$_{e}$", "\n\nGLU$_{i}$", "\n\nGLU$_{a}$", "\n\nGLN$_{e}$", - "\n\nGLN$_{i}$","\n\nGLN$_{a}$", "\n\nGABA$_{i}^{c}$", "\n\nGABA$_{i}^{v}$"] -# plt.title("Concentration control") -plt.xticks(np.arange(len(xticks)), xticks, rotation=60, ha="left") -plt.yticks(np.arange(len(yticks)), yticks, rotation=0, va="center") - -plt.xlabel("Parameter") -# plt.ylabel("Metabolite pool") -plt.ylabel("Metabolite pool") -plt.grid(alpha=0.5) - -for spine in plt.gca().spines.values(): - spine.set_visible(True) # make sure they're visible - spine.set_linewidth(0.5) - spine.set_color("black") - -# Save -plt.tight_layout(rect=[0, -0.01, 1, 1.05]) -plt.savefig(OUTDIR + f"fig_mca_concentration.pdf", dpi=300, transparent=True) diff --git a/reproduce_paper/results/model_fluxes.csv b/reproduce_paper/results/model_fluxes.csv new file mode 100644 index 0000000..7cfff6e --- /dev/null +++ b/reproduce_paper/results/model_fluxes.csv @@ -0,0 +1,56 @@ +Ketone transport,Ketone utilization,PMAS in e,GLU synaptic flux,PMAS in i,GAD67,GAD65,GT,GABA synaptic flux,Anaplerosis,GDH,GS,GLN uptake in e,GLN uptake in i,GLN efflux,CMRglc,CMRglc(ox)N,VPMAS,timestamp +0.00047525787095731585,0.0016058375876988565,0.1903924673331909,0.19120749129283565,0.04775344027329294,0.0480635075000875,0.047930456183443264,0.04808642047901576,0.047933188760640996,0.06,0.04830859634813485,0.25082826844297673,0.1903924673331909,0.04775344027329294,0.011692992251260028,0.4088418438140272,0.27903455840307356,0.23814590760648385,2.3280022923645447 +-2.6806125025118402e-5,0.006455929560399523,0.18789429852008926,0.1911022576317843,0.04720890740057973,0.047891554769279636,0.047905933888784204,0.04808346234352011,0.04792162210434863,0.06,0.048312602937517324,0.25068906164441623,0.18789429852008926,0.04720890740057973,0.011705410187984155,0.40526090365561157,0.27659056674495486,0.235103205920669,5.179572512606683 +-0.001318629806831637,0.018983936665619115,0.18144141539607134,0.19044021820777693,0.04580235434204504,0.04705420442237999,0.0477843292361748,0.048048131732673596,0.047832739233350186,0.06,0.04834804853647534,0.24983972580126562,0.18144141539607134,0.04580235434204504,0.011782739264149217,0.3960111731349309,0.27027762566459035,0.22724376973811639,10.850951381606137 +-0.0024616505084042156,0.030127596882595253,0.17570157694052352,0.1892626646400937,0.04455122563260099,0.045945985914740625,0.04761758233281625,0.04795558672817498,0.04768230966912087,0.06,0.04842342615315086,0.24837533934896028,0.17570157694052352,0.04455122563260099,0.011917900173324385,0.38778353938956733,0.2646622656333797,0.22025280257312452,16.257114187361367 +-0.0035638712911610455,0.04092611532541566,0.1701395129315527,0.18724081659138864,0.04333884691858633,0.044569173180501756,0.04740062731178499,0.047745684169660574,0.0474657207009989,0.06,0.048565716192344704,0.2459431493921732,0.1701395129315527,0.04333884691858633,0.012144099979339944,0.37981073216014977,0.2592208246993022,0.21347835985013902,23.072412575098497 +-0.0043500165280622154,0.04865980382353389,0.16615607106591657,0.1847136999496438,0.04247056490697648,0.04334527947090573,0.047197998804418335,0.04741259802169342,0.04725302589107104,0.06,0.04874869314401736,0.24300391696814055,0.16615607106591657,0.04247056490697648,0.012417692567648917,0.3741007629736957,0.2553237707295473,0.20862663597289305,30.421922748109864 +-0.004851411472858034,0.053606158933074535,0.16360831936497308,0.1818239665227219,0.04191522431293862,0.0423083225713617,0.04701860313622252,0.046934596877139294,0.04706185403089869,0.06,0.04894983457760954,0.23974006319514848,0.16360831936497308,0.04191522431293862,0.01271944123584907,0.37044874941411043,0.25283127147513035,0.2055235436779117,38.590793844286964 +-0.005086816965364112,0.055932223121792954,0.1624102181647822,0.17890561657831724,0.04165407082926707,0.041464343884338596,0.04686703996765089,0.04632010735365321,0.0469011374863776,0.06,0.049132980884190464,0.23651876898710203,0.1624102181647822,0.04165407082926707,0.013013280280199276,0.36873135998906076,0.25165915319253396,0.20406428899404927,47.35223890575253 +-0.005107504099260205,0.05613674924524258,0.1623048715474516,0.1761100154797131,0.04163110813120407,0.04073038220773056,0.04673096509717145,0.04554699975588004,0.04675896507433156,0.06,0.049281347767310225,0.23348220611697212,0.1623048715474516,0.04163110813120407,0.01328529977879615,0.3685803534075303,0.25155609120063943,0.20393597967865568,56.98591403314276 +-0.0049582928107510076,0.05466196667037416,0.16306449752060412,0.17373848207408063,0.041796685938107675,0.04007916470337468,0.046606751809059,0.04464751180116694,0.046631432199484,0.06,0.04937975203608478,0.23093010351312593,0.16306449752060412,0.041796685938107675,0.013509171814657961,0.3696692210366399,0.25229924335750675,0.2048611834587118,67.12458142097785 +-0.0046860671089818794,0.05197380716258202,0.16444910554140202,0.1718898457716877,0.042098492837302263,0.03946102096523707,0.0464856930563757,0.043630046287109876,0.046508686480094465,0.06,0.049431037558765094,0.2289450086361457,0.16444910554140202,0.042098492837302263,0.01367919674953116,0.3716539541897144,0.2536538237344801,0.2065475983787043,77.83192277720445 +-0.004336037593535644,0.048522054308434065,0.16622702261582104,0.17064330358137458,0.04248603041558835,0.038862331587414964,0.04636539925484909,0.04253933931234809,0.04638795452481175,0.06,0.04944035531692637,0.22759625002447342,0.16622702261582104,0.04248603041558835,0.013790957011373495,0.374202466770902,0.2553931835711406,0.20871305303140938,88.90236214333345 +-0.003943658545667363,0.044658950594139844,0.16821681688676413,0.16998175985067698,0.042919751461493044,0.03828155593034246,0.04624573039018161,0.04142091147028466,0.04626619513441393,0.06,0.04941673980401272,0.22686285869616926,0.16821681688676413,0.042919751461493044,0.013847763290700973,0.3770546896368816,0.2573398256771717,0.21113656834825717,100.25709671195546 +-0.00354269487094243,0.040718159549138086,0.1702466260827517,0.16984269603235944,0.04336219467299823,0.03774264980460802,0.04613197018922855,0.04034589318387074,0.0461577973508666,0.06,0.049360381721600156,0.22666903320239237,0.1702466260827517,0.04336219467299823,0.013856000137878223,0.3799642709372778,0.2593256149146921,0.21360882075574994,111.62881659645609 +-0.0031880989763214203,0.037238824060385194,0.17203875032642457,0.1700799535233795,0.0437528290278698,0.037308941430254136,0.046038450995825986,0.03945151665800413,0.04604708769200504,0.06,0.0493074208156457,0.22690798835276965,0.17203875032642457,0.0437528290278698,0.013829290786980862,0.3825331484444271,0.26107887381332145,0.21579157935429438,121.88040170262502 +-0.0028760927641743738,0.03418182831531007,0.17361333725696454,0.17053522096095877,0.04409604616360814,0.036969592877851336,0.04596402019889286,0.038719508240246114,0.045986720618530444,0.06,0.04922364772696311,0.2273458141440257,0.17361333725696454,0.04409604616360814,0.013782095634452854,0.3847902023084077,0.26261931307548825,0.21770938342057267,131.25997424234083 +-0.002588896915125244,0.03137159136890959,0.1750608245080185,0.171144812380457,0.0444115590252284,0.03670327990716568,0.0459048170786221,0.03810497804646762,0.04591376055406076,0.06,0.049153587914392394,0.22799263294859373,0.1750608245080185,0.0444115590252284,0.013719556599328948,0.38686506818572947,0.26403540903676037,0.21947238353324688,140.33753265450468 +-0.0023148871206765853,0.02869365449057545,0.1764401670917395,0.17189064678696164,0.04471221820231446,0.036502323385781424,0.045859673966302625,0.03758615530952853,0.045871085771511035,0.06,0.049062746943284274,0.22877427266895753,0.1764401670917395,0.04471221820231446,0.013643012748402534,0.38884225372433473,0.2653848381668584,0.22115238529405398,149.5266464142187 +-0.002049868219474749,0.026106622270965864,0.17777268681423022,0.17276218008277403,0.045002671268828426,0.03636985804886508,0.04582969313991462,0.03716202946362792,0.04583327798732971,0.06,0.04897061630111863,0.22971695890968472,0.17777268681423022,0.045002671268828426,0.013552892405068654,0.39075232215570854,0.2666884598712711,0.22277535808305865,159.0424034776588 +-0.0017992544845451266,0.023662949899139078,0.17903136523581412,0.17372203309083656,0.04527702893958818,0.03631283739845533,0.045816732483676514,0.036845667505089136,0.04582036247630334,0.06,0.048866239729120575,0.2307592652147786,0.17903136523581412,0.04527702893958818,0.013452681181551487,0.3925565445456219,0.2679198416523869,0.2243083941754023,168.7569889777078 +-0.001566311017852597,0.02139396137529221,0.18020006809742267,0.17473370533907384,0.045531774384828805,0.036330459595900214,0.04582074152715086,0.036638374583920344,0.04581667832704406,0.06,0.0487644316312751,0.23188340949614616,0.18020006809742267,0.045531774384828805,0.013345824600345734,0.3942317936274456,0.2690631991507316,0.22573184248225148,178.56801146610815 +-0.0013532309981654864,0.019320462607231886,0.18126807876310258,0.1757618553944741,0.045764571669312144,0.03641561819614148,0.04584007023742147,0.03653278244766552,0.0458375797874498,0.06,0.04865493413261996,0.23302735352013582,0.18126807876310258,0.045764571669312144,0.013235999530810328,0.3957627078964462,0.2701080481393245,0.22703265043241472,188.3642554633061 +-0.0011599710874030734,0.017441493129331326,0.18223589194620954,0.1767814332037839,0.04597552862709782,0.03655764664860737,0.0458721424771434,0.03651518474145624,0.045861965567419456,0.06,0.04855425567851022,0.23418705320783229,0.18223589194620954,0.04597552862709782,0.013125819373377909,0.3971499964876159,0.2710548726027978,0.22821142057330734,198.09287364881675 +-0.000984809423861322,0.015739844824217288,0.18311237114508458,0.17777923402356072,0.046166577259813824,0.03674595302925125,0.04591435101592075,0.036570540621226326,0.04590730948857384,0.06,0.04844767546831138,0.23531823202669838,0.18311237114508458,0.046166577259813824,0.013016912727640252,0.3984063645700318,0.2719123438190467,0.2292789484048984,207.76919376142177 +-0.0008259046845793044,0.014197246393597113,0.1839069274776667,0.17874708603256728,0.046339768935269976,0.036971087186205655,0.04596435041482389,0.036685261366270386,0.04595063116629142,0.06,0.0483522647284434,0.23643831818667327,0.1839069274776667,0.046339768935269976,0.012910234594757328,0.3995453022849471,0.2726896688094764,0.23024669641293669,217.4207924223956 +-0.0006815652333935964,0.01279696437903948,0.18462817995658973,0.1796797459220473,0.04649698236385804,0.0372249337904237,0.04602012928546877,0.036847666258120394,0.04601026367906309,0.06,0.048252368289004026,0.23751087493214265,0.18462817995658973,0.04649698236385804,0.012806640675989873,0.4005791643492668,0.27339527966837457,0.23112516232044777,227.07720359071573 +-0.0005505935009987422,0.011527124839285634,0.18528224457097478,0.18057178981237265,0.046639550664904494,0.03749999321328285,0.04607986664268619,0.037047425986882127,0.04606438473605368,0.06,0.04816476427316267,0.2385562299510044,0.18528224457097478,0.046639550664904494,0.012706828854564737,0.40151671758088453,0.2740351597489537,0.23192179523587927,236.7443228690372 +-0.0004318620564750923,0.010376584067623627,0.18587486118322577,0.1814195697158347,0.046768724972338895,0.03778983325622383,0.046142037342177554,0.037275664582341134,0.04613098743064714,0.06,0.04807423441797439,0.23954057891424463,0.18587486118322577,0.046768724972338895,0.012611497171029358,0.40236618964787396,0.27461492443467395,0.23264358615556466,246.42737933902572 +-0.00032451402387394793,0.009336861123490088,0.18641039812725085,0.1822193715529804,0.04688545746417238,0.03808847894461068,0.04620527875421175,0.037524405312063,0.046189202236679036,0.06,0.04799730402163022,0.24048635846453253,0.18641039812725085,0.04688545746417238,0.012521150373513932,0.40313384223368687,0.2751388473244913,0.23329585559142324,256.1123923712603 +-0.00022753201158790635,0.008397954255796564,0.1868940070708505,0.18297025333850248,0.04699087106390472,0.03839141699058318,0.04626859709527225,0.03778731919524017,0.04625745938293606,0.06,0.0479183664247861,0.24136288384513654,0.1868940070708505,0.04699087106390472,0.012436165153799294,0.40382705987295237,0.27561196836329,0.23388487813475523,265.8016211130471 +-0.00014008916942835766,0.007551736492821401,0.18732987402687032,0.18367128045770328,0.047085878208512996,0.03869460581536283,0.04633114511147525,0.03805887191764196,0.046315280675610616,0.06,0.04785393886633154,0.2421969203218732,0.18732987402687032,0.047085878208512996,0.012356692968628746,0.4044518429053688,0.2760383827829142,0.2344157522353833,275.4836298929331 +-6.122837554452795e-5,0.006788845236163164,0.18772282145052716,0.18432384691701462,0.04717153006210549,0.03899542388879643,0.046392406288928524,0.03833511092950651,0.046381780714295204,0.06,0.04778774881136155,0.24296015074976368,0.18772282145052716,0.04717153006210549,0.01228280552005208,0.4050151039636844,0.2764228084552146,0.23489435151263266,285.1672855488159 +9.79627425539343e-6,0.006101982940657119,0.18807660813646776,0.1849287806030658,0.047248645939708546,0.0392914084195991,0.046451922313605426,0.03861249763215367,0.04643675324288015,0.06,0.047736322058832846,0.24368241299969146,0.18807660813646776,0.047248645939708546,0.01221440144061313,0.40552223100191864,0.2767689226588095,0.2353252540761763,294.8451494681176 +7.379562313319173e-5,0.005483241496848935,0.18839530736670815,0.185488442532759,0.04731811370604368,0.03958117884835742,0.04650947139423093,0.03888871177025467,0.04649967053792901,0.06,0.04768297103410179,0.24433635581329172,0.18839530736670815,0.04731811370604368,0.012151398025575813,0.4059790627684665,0.2770807103394784,0.23571342107275184,304.5270978868357 +0.00013139365976945477,0.004926534034660728,0.18868205434550517,0.18600443308844816,0.04738061675038853,0.039863280163218766,0.04656482785458491,0.03916153356794987,0.046550578888710066,0.06,0.0476441791993631,0.24495337740839654,0.18868205434550517,0.04738061675038853,0.012093557523258712,0.4063900933466619,0.27736123870909674,0.23606267109589368,314.2057002718702 +0.0001832506706484142,0.004425434522286961,0.18894015897138003,0.18647939158612375,0.047436876540620114,0.04013706684272493,0.046617932218609215,0.03942965421022202,0.04660910034105258,0.06,0.04760300467571873,0.2455063602808034,0.18894015897138003,0.047436876540620114,0.012040727717450612,0.4067600672311411,0.2776137458852538,0.23637703551200015,323.888603206904 +0.00022988702209508015,0.003974878793407377,0.18917222967782987,0.18691530780951307,0.04748746164431406,0.04040175662529613,0.04666870086310136,0.03969173966448902,0.0466554429486018,0.06,0.0475759161178189,0.2460276610118097,0.18917222967782987,0.04748746164431406,0.011992587873756648,0.40709272341832897,0.2778407837330095,0.23665969132214393,333.5685026590654 +0.0002718420268887561,0.0035696274276773264,0.18938096516759628,0.18731486842500392,0.04753296030528668,0.040657175109383305,0.04671716747290383,0.03994714029770808,0.046709336427110004,0.06,0.04754576530170995,0.24649006032459217,0.18938096516759628,0.04753296030528668,0.011948952876451818,0.4073919302989332,0.2780449924290219,0.23691392547288295,343.2515331952727 +0.0003095518583297752,0.0032054437649170677,0.1895685476451389,0.18768023623273164,0.04757384818477314,0.04090300335117972,0.046763336095660114,0.04019513351162704,0.04675106063184903,0.06,0.04752906226655263,0.24692654043544596,0.1895685476451389,0.04757384818477314,0.011909454953210535,0.4076608159015893,0.27822850685283473,0.23714239582991203,352.931661562523 +0.0003434566817597926,0.0028780577827392136,0.18973717649900945,0.18801399152934398,0.04761060468967341,0.04113935860458375,0.046807289961887315,0.0404354923193943,0.04680042503475399,0.06,0.047508506888403655,0.2473095668443753,0.18973717649900945,0.04761060468967341,0.011873903900314184,0.4079025328870572,0.27839347869541653,0.23734778118868285,362.6141868395043 +0.00037391900697976866,0.0025839536135261265,0.18988866267029586,0.1883182811598535,0.04764362455590274,0.04136621456661214,0.04684908137253139,0.04066789266673981,0.04683773684698891,0.06,0.047500680602408585,0.2476724238796602,0.18988866267029586,0.04764362455590274,0.011841916855724315,0.4081196771053979,0.27854167962443405,0.2375322872261986,372.2943370521789 +0.0004012959054878074,0.0023196713305575004,0.19002478828623953,0.1885954715770807,0.047673296238741364,0.04158384750337654,0.04688881425345634,0.04089235922571501,0.04688284918317143,0.06,0.047488194545823134,0.24798715094422197,0.19002478828623953,0.047673296238741364,0.011813319076241923,0.40831480310415563,0.2786748531185862,0.2376980845249809,381.9765531697594 +0.00042588582640948396,0.002082319750208187,0.19014704252783055,0.18884757290033266,0.04769994433915227,0.041792400533258126,0.0469265634811906,0.041108813859085824,0.0469160762997886,0.06,0.047487712052331495,0.24828709878592237,0.19014704252783055,0.04769994433915227,0.011787734484425433,0.4084900455140976,0.2787944560633716,0.2378469868669828,391.65690525132413 +0.0004479775323769195,0.001869103447580355,0.19025686525298044,0.18907669262265717,0.04772388270700854,0.041992226431405315,0.046962437537825585,0.04131742397118359,0.0469572937777893,0.06,0.04748180133093513,0.24854384928570306,0.19025686525298044,0.04772388270700854,0.011765014927532354,0.4086474682650076,0.27890189709086766,0.23798074795998897,401.33900773681853 +0.0004678154591330237,0.0016776566707203853,0.19035547500512792,0.18928464776769252,0.04774537695164596,0.04218356084286929,0.04699651936143661,0.04151825986980625,0.046986808194068704,0.06,0.04748721234224424,0.2487906888546775,0.19035547500512792,0.04774537695164596,0.011744803230565443,0.40878881804874695,0.2789983683182698,0.23810085195677388,411.01955731786074 +0.00048563302766502995,0.0015057214490399693,0.19044403481205388,0.18947328053153478,0.047764680581496514,0.04236678619075904,0.04702891398924714,0.04171156676854562,0.04702451022202175,0.06,0.047486498439349156,0.24899875228452686,0.19044403481205388,0.047764680581496514,0.011726982526796471,0.40891576197958923,0.2790850075510697,0.2382087153935504,420.70155986838836 +0.0005016297107930475,0.0013513686163809649,0.1905235383420762,0.18964419415912168,0.047782010189234336,0.04254218254862594,0.047059704364426774,0.041897500875609214,0.04705068857166098,0.06,0.04749649553229368,0.249201192383045,0.1905235383420762,0.047782010189234336,0.011711220550105388,0.4090297244075483,0.2791627869081517,0.23830554853131053,430.3822141217006 +0.0005159941260225197,0.0012127747204613288,0.19059492481277027,0.18979897614198935,0.047797570498752266,0.042710134069669885,0.047088988158629545,0.04207634473506142,0.047085245673922034,0.06,0.04749976032100973,0.24936872905750188,0.19059492481277027,0.047797570498752266,0.011697432815072564,0.4091320516315275,0.2792326252385175,0.23839249531152254,440.06410351469276 +0.0005288885583984937,0.0010883711670193934,0.1906590021703126,0.189939015747112,0.04781153762034764,0.042870935064123115,0.04711684397838746,0.042248299902677144,0.047108447613478044,0.06,0.04751320800933556,0.2495343465157254,0.1906590021703126,0.04781153762034764,0.011685314068901715,0.4092239017825921,0.2792953129666191,0.23847053979066024,449.7448208888509 +0.0005404652961790222,0.0009766864428990061,0.19071652835667516,0.19006566487890178,0.04782407676476093,0.04302495448252195,0.04714336036758931,0.04241366031531509,0.04714020610527426,0.06,0.047519410364025,0.24966839373263258,0.19071652835667516,0.04782407676476093,0.011674811079771339,0.40930636131448045,0.2793515915971329,0.2385406051214361,459.4266189596274 +0.0005508560442727726,0.0008764481251350169,0.19076815876713443,0.19018011082865302,0.04783533079034131,0.0431724824924384,0.04716860948517959,0.04257264690528703,0.04716076302242269,0.06,0.047535351942104294,0.24980367397353317,0.19076815876713443,0.04783533079034131,0.011665645209713471,0.4093803696878295,0.2794021023119436,0.23860348955747573,469.10738380451403 +0.0005601836913123134,0.0007864692533460497,0.1908145047771908,0.1902834942376565,0.047845432960323765,0.043313863977169687,0.047192670379575634,0.04272554959142599,0.04719003816734677,0.06,0.04754362596374097,0.24991021385646847,0.1908145047771908,0.047845432960323765,0.011657791950260125,0.40944680326418464,0.279447443227806,0.23865993773751457,478.78911556306593 +0.0005685549213729635,0.0007057195194862127,0.1908560970743948,0.1903768189535973,0.04785449895016417,0.043449374467551864,0.04721560798893902,0.04287259254986206,0.04720824894209255,0.06,0.04756127555276009,0.2500206383579105,0.1908560970743948,0.04785449895016417,0.0116509974140222,0.4095064227448848,0.27948813352338386,0.23871059602455896,488.4699123056489 +0.0005760688632734977,0.0006332417631528852,0.19089342866965825,0.19046104026511917,0.04786263622283573,0.04357933089748263,0.047237492238520834,0.04301405296479589,0.04723532298445146,0.06,0.0475709179515754,0.25010469947581115,0.19089342866965825,0.04786263622283573,0.011645262619425993,0.4095599348245952,0.2795246555177862,0.23875606489249399,498.15158896429966 +0.0005774132918350931,0.0006202739626317107,0.19090010808007332,0.19047615619229138,0.047864092152684126,0.043603532496091515,0.047241555533928,0.04304044450540069,0.04723796752087256,0.06,0.0475749052482255,0.2501234350675255,0.19090010808007332,0.047864092152684126,0.011644257079321024,0.4095695092652272,0.27953119007351757,0.23876420023275746,500.0 diff --git a/reproduce_paper/results/model_mca_concentration.csv b/reproduce_paper/results/model_mca_concentration.csv new file mode 100644 index 0000000..12abff6 --- /dev/null +++ b/reproduce_paper/results/model_mca_concentration.csv @@ -0,0 +1,10 @@ +GLN_p,c_vpmas,Tmax_gln,KI_GS,KT_ket,Vmax_GT,KM_GAD65,V_anapl,V_ATP,V_out0,KM_GS,k_ket,KM_GDH,Vmax_GS,GLN_ratio,Vmax_GAD67,k_syn_i,Tmax_ket,ATPglc,KM_GT,KM_GAD67,k_syn_e,Vmax_GAD65,KT_gln,h_GDH,ATPket,Vmax_GDH,state +NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,KET_b(t) +0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,GLN_e(t) +0.0,0.999999999999955,0.0,0.0,0.0,0.0,0.0,0.0,1.4959093516347277,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,-1.4960589575304735,0.0,0.0,-1.0001000100009545,0.0,0.0,0.0,0.0,0.0,GLU_e(t) +0.06081445392049228,0.17278439832680503,0.002605986935952624,-0.060296769516393206,0.0,0.001667701587154338,0.0009882897034915144,0.5500118443544302,0.2487379569568663,-0.10138563865477676,0.12011612376192589,0.0,0.8799534969642383,-0.1799255275808817,-0.06028488385847596,-0.0016586704883029659,-0.00017294257769141116,0.0,-0.24858988443011057,-0.0006628810452739565,0.0009586089729876369,5.185853500141592e-5,0.001660763042142853,-0.0005301281933141901,-0.002827768615007273,0.0,-0.44288440037089377,GLU_a(t) +0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,GLN_i(t) +-3.184376928761778e-14,14.297283779206886,-2.122917952507852e-14,2.122917952507852e-14,0.0,6.07216696639272e-6,0.9999979268866238,2.122917952507852e-14,18.722035894932905,-1.061458976253926e-14,-1.061458976253926e-14,0.0,-2.122917952507852e-14,1.061458976253926e-14,0.0,-5.893486577448119e-6,2.122917952507852e-14,0.0,-18.780985355920617,-5.529781777695539e-7,1.995129203565155e-6,-1.061458976253926e-14,-14.334330966559698,1.061458976253926e-14,2.122917952507852e-14,0.0,0.0,GLU_i(t) +-2.9878346779710437e-14,47.649892624932534,-2.9878346779710437e-14,-7.469586694927609e-15,0.0,-5.773816619124541,2.870159980613984,2.9878346779710437e-14,70.01234159893866,-7.469586694927609e-15,-1.4939173389855218e-14,0.0,-1.4939173389855218e-14,-7.469586694927609e-15,0.0,5.767735591040449,0.0,0.0,-70.91888596361385,0.9999981976381839,-2.8715213046638675,1.4939173389855218e-14,-48.08669190975708,-2.240876008478283e-14,-2.240876008478283e-14,0.0,-1.4939173389855218e-14,GABA_i(t) +-3.113686741893949e-6,0.9948970484544969,-1.617758767585503e-7,3.150552790020523e-5,0.0,-0.003942517075640261,-0.0023374625008876454,-1.7618662551555943e-5,1.29068248269248,5.10132595208432e-6,-0.00011109735666857862,0.0,-4.360245180954103e-5,-0.00032661161932168425,3.0895631267810253e-6,0.003921189759779586,-0.9996820346196911,0.0,-1.2903670445945512,0.0015674167163353377,-0.0022667026482348744,-0.00012262467135280733,-0.001245121193565172,3.3315936969961176e-8,-9.203046448782115e-7,0.0,1.4693887172371982e-5,GABA_v(t) +0.19063389198417954,-2.3107472025881735,0.008168912028112709,0.8110357741487492,0.0,-0.0001594176012200251,-9.488435329662688e-5,1.0096997783781823,-3.3703181353213534,-0.3178121958229346,-1.6155414875733638,0.0,1.6153901993550392,2.422325709347045,-0.18897446132743367,0.00015856427695928803,1.6483288450382186e-5,0.0,3.370170951067615,6.34803987455e-5,-9.174903174635032e-5,-4.967698545114672e-6,0.0009993538244214184,-0.0016617770558616039,-0.005190593862305346,0.0,-0.8130352414338525,GLN_a(t) diff --git a/reproduce_paper/results/model_multirun.h5 b/reproduce_paper/results/model_multirun.h5 new file mode 100644 index 0000000..67d45c1 Binary files /dev/null and b/reproduce_paper/results/model_multirun.h5 differ diff --git a/reproduce_paper/results/model_timeseries.csv b/reproduce_paper/results/model_timeseries.csv new file mode 100644 index 0000000..a4e0d56 --- /dev/null +++ b/reproduce_paper/results/model_timeseries.csv @@ -0,0 +1,56 @@ +timestamp,KET_b,GLN_e,GLU_e,GLU_a,GLN_i,GLU_i,GABA_i,GABA_v,GLN_a +2.3280022923645447,0.009912577701844793,0.28,10.010863418473072,0.8025758095337949,0.07,0.13061693381935985,0.7431500889745046,0.10000665295355934,3.8323004212629437 +5.179572512606683,0.03985141703950323,0.28,10.00535380271122,0.8026093073824924,0.07,0.12968347735665536,0.7428870844441885,0.09998252055987614,3.838974355236782 +10.850951381606137,0.11718479423221675,0.28,9.97069205276319,0.802905719564575,0.07,0.12523179329096823,0.7397577265642256,0.09979707747412933,3.880739075429039 +16.257114187361367,0.18597282026293366,0.28,9.90904003351276,0.803536436336686,0.07,0.11956953366594039,0.7316630056528797,0.09948322484690354,3.95459607708033 +23.072412575098497,0.2526303415149115,0.28,9.80318411473239,0.8047284389145366,0.07,0.11287549553486549,0.7138344160275529,0.0990313388295408,4.08070160516872 +30.421922748109864,0.3003691594045302,0.28,9.670874342913288,0.8062640119665416,0.07,0.10721941386974505,0.6869669903553556,0.09858757749023793,4.237564136138565 +38.590793844286964,0.33090221563626254,0.28,9.519579399095388,0.8079556186338602,0.07,0.1026294043809476,0.6511864166680629,0.09818872111599976,4.416329845387623 +47.35223890575253,0.34526063655427747,0.28,9.366786208288861,0.8094992157992336,0.07,0.09902268273410567,0.6094190406343726,0.09785340598034133,4.596513215566945 +56.98591403314276,0.34652314348915175,0.28,9.220419658623724,0.8107520486873591,0.07,0.09597582689338333,0.5626330454828172,0.09755678087696967,4.7689637036778265 +67.12458142097785,0.3374195473479886,0.28,9.096255605972809,0.8115841716817188,0.07,0.093339545776265,0.5149169652199562,0.09729069935214688,4.915158291443117 +77.83192277720445,0.3208259701393952,0.28,8.999468365009827,0.8120182276551146,0.07,0.09089343594116218,0.46804378167942495,0.0970346056334122,5.028868937246041 +88.90236214333345,0.2995188537557658,0.28,8.934204375988198,0.8120971165336468,0.07,0.08857485319735171,0.4245605512519989,0.09678271338370906,5.104912247212971 +100.25709671195546,0.2756725345317274,0.28,8.899568578569475,0.8118971922693572,0.07,0.08637168794298021,0.38584267860200205,0.09652867751807621,5.1439674417582015 +111.62881659645609,0.2513466638835684,0.28,8.89228775038531,0.8114202980305888,0.07,0.0843667130313617,0.35319958425151493,0.09630251898782934,5.149653217126757 +121.88040170262502,0.22986928432336537,0.28,8.904709608553901,0.8109724344775292,0.07,0.08277990711090298,0.328894446615521,0.09607153701649288,5.131237203565623 +131.25997424234083,0.21099894021796337,0.28,8.928545600050198,0.810264567121209,0.07,0.08155462194854494,0.3106625444873742,0.09594558860532118,5.098844582064923 +140.33753265450468,0.19365179857351597,0.28,8.96046138117576,0.8096730974304414,0.07,0.08060288442415095,0.29638226235507814,0.09579336648041052,5.056209930534245 +149.5266464142187,0.17712132401589784,0.28,8.999510302982285,0.8089068900257979,0.07,0.07989036458398538,0.28498885885148983,0.09570433083978935,5.00447149779497 +159.0424034776588,0.16115198932694977,0.28,9.045140318469844,0.8081306108385282,0.07,0.07942332054631226,0.2760924679081085,0.09562544958758545,4.9441735552014165 +168.7569889777078,0.1460675919699943,0.28,9.0953944026616,0.807252121707147,0.07,0.07922291741538338,0.26968662365905854,0.09559850297580501,4.877895203138227 +178.56801146610815,0.13206148997093956,0.28,9.148361536077164,0.8063962366535999,0.07,0.07928481091893697,0.26559119903808986,0.09559081645533915,4.808100412088654 +188.3642554633061,0.11926211485945609,0.28,9.202191381909639,0.8054767850304839,0.07,0.07958442572553244,0.263535171637897,0.09563442475996203,4.737293901346492 +198.09287364881675,0.10766353783537855,0.28,9.255572419046278,0.804632362807967,0.07,0.08008603966508226,0.2631944661295162,0.09568530266517725,4.667184616716573 +207.76919376142177,0.09715953595195856,0.28,9.307813299662865,0.8037394496153747,0.07,0.080754811182485,0.26426807064383284,0.0957799071324303,4.598779663162477 +217.4207924223956,0.08763732341726613,0.28,9.358486179715564,0.8029409847640605,0.07,0.08155998651109797,0.2665106208893824,0.0958702924395815,4.532619713128474 +227.07720359071573,0.07899360727802147,0.28,9.407316540421325,0.802105855240941,0.07,0.08247526087947171,0.26972649949022576,0.09599470828095784,4.469158272556672 +236.7443228690372,0.0711550916005286,0.28,9.454020409024746,0.801374219952779,0.07,0.08347600258901976,0.27374978941847267,0.09610762515346064,4.408732931988193 +246.42737933902572,0.06405298807175078,0.28,9.498406791404959,0.800618859886847,0.07,0.08454077141045659,0.27844068917491027,0.09624658341466126,4.351667522525605 +256.1123923712603,0.05763494520672893,0.28,9.540281233140336,0.799977533829023,0.07,0.08564905701158333,0.2836708998603734,0.09636804138677037,4.298160568144153 +265.8016211130471,0.05183922380121335,0.28,9.579594415628403,0.7993200087071379,0.07,0.08678503186139322,0.2893376938926532,0.09651045145615701,4.248331390962516 +275.4836298929331,0.046615657363095066,0.28,9.616297406162476,0.7987837456141139,0.07,0.08793399342309023,0.2953460032826026,0.09663108841145548,4.202168996999551 +285.1672855488159,0.04190645207508126,0.28,9.650463189372493,0.798233182471211,0.07,0.08908607096267349,0.30162654405061917,0.09676983249383518,4.159622066784984 +294.8451494681176,0.03766656136208098,0.28,9.682135110108158,0.7978056766446875,0.07,0.09023158767858024,0.30811143077703124,0.09688452585620728,4.1205477790906135 +304.5270978868357,0.0338471697336354,0.28,9.711436781819843,0.79736241116133,0.07,0.0913647186689452,0.3147540224855815,0.09701579498837681,4.084823464489126 +314.2057002718702,0.030410703917658814,0.28,9.73845199415959,0.7970402609220215,0.07,0.09247911791255388,0.32150424082300627,0.0971220089478616,4.05224796794798 +323.888603206904,0.027317497051154077,0.28,9.763318931210668,0.7966984615921823,0.07,0.09357145318804849,0.32832932446102786,0.09724410669946293,4.022677821690017 +333.5685026590654,0.024536288848193685,0.28,9.786141770131575,0.7964736704794404,0.07,0.09463774061869618,0.3351918979471366,0.09734079480200668,3.9958838470562075 +343.2515331952727,0.02203473720788473,0.28,9.807061174083975,0.7962235395319368,0.07,0.09567636656317316,0.3420690308460419,0.09745323686023369,3.971720685676129 +352.931661562523,0.01978668990689548,0.28,9.82619037867705,0.7960850038226956,0.07,0.09668510572690381,0.34893340486164504,0.09754028923815779,3.949948966957447 +362.6141868395043,0.017765788782340825,0.28,9.843664477976125,0.7959145483928866,0.07,0.09766351310566557,0.35576935556853195,0.09764328194190276,3.9304339450743604 +372.2943370521789,0.01595033094769214,0.28,9.85959587224364,0.7958496581372023,0.07,0.09861057417751128,0.362557193096389,0.09772112841015837,3.912940620978751 +381.9765531697594,0.014318958830601854,0.28,9.874108459533021,0.7957461427198277,0.07,0.09952656672018523,0.369286364925915,0.09781524970409228,3.8973529449707303 +391.65690525132413,0.012853825618569053,0.28,9.88730748169281,0.7957421428753813,0.07,0.10041125688200317,0.3759427551253289,0.0978845739615869,3.8834491428905187 +401.33900773681853,0.01153767560234787,0.28,9.89930327867315,0.7956931448826915,0.07,0.10126534773572916,0.3825194304277472,0.09797056911702337,3.8711350343672932 +411.01955731786074,0.010355905374817193,0.28,9.910190982601703,0.7957380003265958,0.07,0.10208909964287981,0.3890064206134681,0.09803214728576821,3.86020596071479 +420.70155986838836,0.009294576845925736,0.28,9.920067043535854,0.795732082175752,0.07,0.10288345893868103,0.3953994775856955,0.09811080789071927,3.8505898253510153 +430.3822141217006,0.008341781582598549,0.28,9.929015401001136,0.7958149604549726,0.07,0.10364898707413217,0.40169174903129684,0.0981654257702086,3.8421002064759024 +440.06410351469276,0.007486263706551413,0.28,9.937119169737663,0.795842028121212,0.07,0.10438675360320893,0.40788100303725733,0.09823752487778434,3.8346859312847617 +449.7448208888509,0.006718340537156749,0.28,9.944451086236231,0.7959535293385223,0.07,0.10509748880307776,0.4139627271388254,0.0982859328468142,3.8281783633734334 +459.4266189596274,0.0060289286598704086,0.28,9.951081930832554,0.7960049610778582,0.07,0.10578230293681705,0.41993617518031173,0.09835219300078084,3.8225453927682596 +469.10738380451403,0.005410173611944549,0.28,9.957073865374504,0.7961371679483963,0.07,0.10644200579350471,0.4257985770826137,0.09839508245863278,3.8176348204717465 +478.78911556306593,0.004854748477444751,0.28,9.962486609301388,0.7962057944741945,0.07,0.10707769657362624,0.4315502796428729,0.09845616141737278,3.8134313770665527 +488.4699123056489,0.0043562933301618065,0.28,9.967372720083628,0.7963522030064272,0.07,0.10769020623571778,0.43718980116778555,0.09849415594010547,3.809797516599346 +498.15158896429966,0.003908899772548674,0.28,9.971782212833466,0.7964322004938831,0.07,0.10828059185504874,0.4427182883498867,0.09855064257135711,3.8067325218450425 +500.0,0.0038288516211833997,0.28,9.972573622633057,0.7964652830859185,0.07,0.10839086261549605,0.4437611299744736,0.09855616006858452,3.8061953023550963 diff --git a/tutorial/NTModelModule.jl b/tutorial/NTModelModule.jl index dfb297b..6a45885 100644 --- a/tutorial/NTModelModule.jl +++ b/tutorial/NTModelModule.jl @@ -27,9 +27,9 @@ mutable struct NTModel :GLN_e => 0.28, # mM :GLN_i => 0.07, # mM :GLU_e => 10.0, # mM - :GLU_i => 0.06, # mM - :GABA_i => 0.6, # mM - :GABA_v => 0.06, # mM + :GLU_i => 0.13, # mM + :GABA_i => 0.73, # mM + :GABA_v => 0.1, # mM :GLN_a => 3.8, # mM :GLU_a => 0.8, # mM :KET_b => 0.0, # mM @@ -64,23 +64,26 @@ mutable struct NTModel ps = Dict([ # Neurotranstmitter cycling :c_vpmas => 1.0, - :k_syn_e => 0.01912, - :k_syn_i => 0.79882, - :Vmax_GAD67 => 0.20769, - :KM_GAD67 => 0.2, - :Vmax_GT => 0.05860, - :KM_GT => 0.13359, - :Vmax_GAD65 => 0.05591, + :k_syn_e => 0.0191, + :k_syn_i => 0.4793, + :Vmax_GAD67 => 0.0959, + :KM_GAD67 => 0.13, + :Vmax_GT => 0.0582, + :KM_GT => 0.1563, + :Vmax_GAD65 => 0.0516, :KM_GAD65 => 0.01, + :Vmax_GDH => 0.096, + :KM_GDH => 0.8, + :h_GDH => 4.0, :Vmax_GS => 0.753, :KM_GS => 0.8, :KI_GS => 3.8, - :h_GS => 4.0, - :Vmax_GDH => 0.0764, - :KM_GDH => 0.8, - :h_GDH => 4.0, - :k_ox => 0.2, - :V_out => 0.012, + :V_anapl => 0.06, + :GLN_p => 0.5, + :GLN_ratio => 1/8, + :V_out0 => 0.012, + :Tmax_gln => 0.018, + :KT_gln => 0.25, # Energy metabolism :V_ATP => 11.94, # mM/min @@ -106,15 +109,18 @@ mutable struct NTModel :KM_GT => "mM", :Vmax_GAD65 => "mM/min", :KM_GAD65 => "mM", - :Vmax_GS => "mM/min", - :KM_GS => "mM", - :KI_GS => "mM", - :h_GS => "-", :Vmax_GDH => "mM/min", :KM_GDH => "mM", :h_GDH => "-", - :k_ox => "-", - :V_out => "mM/min", + :Vmax_GS => "mM/min", + :KM_GS => "mM", + :KI_GS => "mM", + :V_anapl => "mM/min", + :GLN_p => "mM", + :GLN_ratio => "-", + :V_out0 => "mM/min", + :Tmax_gln => "mM/min", + :KT_gln => "mM", # Energy metabolism :V_ATP => "mM/min", @@ -127,7 +133,7 @@ mutable struct NTModel :ATPket => "-", :KET_p => "mM", - ]) + ]); ps_desc = Dict([ # Neurotranstmitter cycling @@ -140,15 +146,18 @@ mutable struct NTModel :KM_GT => "GT Michaelis-Menten constant for GABA", :Vmax_GAD65 => "GAD65 maximum rate", :KM_GAD65 => "GAD65 Michaelis-Menten constant for glutamate", - :Vmax_GS => "GS maximum rate", - :KM_GS => "GS Michaelis-Menten constant for glutamate", - :KI_GS => "GS inhibition constant for glutamine", - :h_GS => "GS Hill exponent for the inhibition by GLN", :Vmax_GDH => "GDH maximum rate", :KM_GDH => "GDH Michaelis–Menten constant for glutamate", :h_GDH => "GDH Hill exponent", - :k_ox => "relative portion of glutamate oxidation", - :V_out => "rate of GLN efflux", + :Vmax_GS => "GS maximum rate", + :KM_GS => "GS Michaelis-Menten constant for glutamate", + :KI_GS => "GS inhibition constant for glutamine", + :V_anapl => "rate of anaplerotic flux", + :GLN_p => "plasma glutamine concentration", + :GLN_ratio => "ratio of extracellular vs. intracellular glutamine", + :V_out0 => "baseline rate of glutamine efflux", + :Tmax_gln => "glutamine transport maximum rate", + :KT_gln => "glutamine transport Michaelis-Menten constant", # Energy metabolism :V_ATP => "rate of ATP synthesis", @@ -160,54 +169,42 @@ mutable struct NTModel :k_ket => "Ketone utilization rate constant", :ATPket => "ATP yield from ketones", :KET_p => "Plasma Ketone concentration", - ]) + ]); # Model definition # ------------- # Custom rate functions - mm_reversible(S1P, S1B, Vmax, Km1) = - Vmax / Km1 * (S1P - S1B) / (1 + S1P/Km1 + S1B/Km1); + mm_reversible(S1P, S1B, Vmax, Km) = Vmax * (S1P - S1B) / (Km + S1P + S1B); CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket) = ((V_ATP - ATPket*KET_b*k_ket)/ATPglc)/0.91 - Vpmas_e(cmrglc, c_vpmas) = c_vpmas*(cmrglc*0.91*0.75*0.83 - 0.077)/0.812; - Vpmas_i(cmrglc, c_vpmas) = c_vpmas*(cmrglc*0.91*0.75*0.17 - 0.011)/0.763; + Vpmas_e(cmrglc) = (cmrglc*0.91*0.75*0.83 - 0.077)/0.812; + Vpmas_i(cmrglc) = (cmrglc*0.91*0.75*0.17 - 0.011)/0.763; # Reaction network model = @reaction_network begin - # Ketone transport + # Ketone metabolism mm_reversible(KET_p, KET_b, Tmax_ket, KT_ket), 0 => KET_b, [description="Ketone transport"] - - # Ketone utilization k_ket*KET_b, KET_b => 0, [description="Ketone utilization"] # Glutamatergic neuron - Vpmas_e(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket), c_vpmas), GLN_e => GLU_e, [description="PMAS in e"] - k_syn_e * GLU_e, GLU_e => GLU_a, [description="GLU cycling"] + Vpmas_e(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket))*c_vpmas, GLN_e => GLU_e, [description="PMAS in e"] + k_syn_e * GLU_e, GLU_e => GLU_a, [description="GLU synaptic flux"] # GABAergic neuron - Vpmas_i(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket), c_vpmas), GLN_i => GLU_i, [description="PMAS in i"] + Vpmas_i(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket))*c_vpmas, GLN_i => GLU_i, [description="PMAS in i"] mm(GLU_i, Vmax_GAD67, KM_GAD67), GLU_i => GABA_i, [description="GAD67"] mm(GLU_i, Vmax_GAD65, KM_GAD65), GLU_i => GABA_v, [description="GAD65"] mm(GABA_i, Vmax_GT, KM_GT), GABA_i => GLU_i, [description="GT"] - # mm(GLU_i, Vmax_GAD65, KM_GAD65), GABA_v => GLU_a, [description="GABA cycling"] - k_syn_i * GABA_v, GABA_v => GLU_a, [description="GABA cycling"] + k_syn_i * GABA_v, GABA_v => GLU_a, [description="GABA synaptic flux"] # Astrocyte - k_ox * k_syn_e * GLU_e + V_out, 0 => GLU_a, [description="Anaplerosis"] + V_anapl, 0 => GLU_a, [description="Anaplerosis"] Vmax_GDH * (GLU_a/KM_GDH)^h_GDH / (1 + (GLU_a/KM_GDH)^h_GDH), GLU_a => 0, [description="GDH"] - Vmax_GS * GLU_a / (KM_GS*(1 + (GLN_a/KI_GS)^h_GS) + GLU_a), GLU_a => GLN_a, [description="GS"] - Vpmas_e(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket), c_vpmas), GLN_a => GLN_e, [description="GLN uptake in e"] - Vpmas_i(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket), c_vpmas), GLN_a => GLN_i, [description="GLN uptake in i"] - V_out, GLN_a => 0, [description="GLN efflux"] - - # Monitoring - CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket), KET_b => KET_b, [description="CMRglc"] - CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket)*0.91*0.75, KET_b => KET_b, [description="CMRglc(ox)N"] - Vpmas_e(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket), c_vpmas) + - Vpmas_i(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket), c_vpmas), - KET_b => KET_b, [description="VPMAS"] - + Vmax_GS * GLU_a / (KM_GS*(1 + GLN_a/KI_GS) + GLU_a), GLU_a => GLN_a, [description="GS"] + Vpmas_e(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket))*c_vpmas, GLN_a => GLN_e, [description="GLN uptake in e"] + Vpmas_i(CMRglc(V_ATP, KET_b, k_ket, ATPglc, ATPket))*c_vpmas, GLN_a => GLN_i, [description="GLN uptake in i"] + V_out0 + mm_reversible(GLN_a*GLN_ratio, GLN_p, Tmax_gln, KT_gln), GLN_a => 0, [description="GLN efflux"] end; diff --git a/tutorial/figures/fig1.png b/tutorial/figures/fig1.png index a984baa..dd4916e 100644 Binary files a/tutorial/figures/fig1.png and b/tutorial/figures/fig1.png differ diff --git a/tutorial/figures/fig2.png b/tutorial/figures/fig2.png index 5d85047..f5a8735 100644 Binary files a/tutorial/figures/fig2.png and b/tutorial/figures/fig2.png differ diff --git a/tutorial/figures/fig3.png b/tutorial/figures/fig3.png index c5ee347..2003262 100644 Binary files a/tutorial/figures/fig3.png and b/tutorial/figures/fig3.png differ diff --git a/tutorial/figures/fig4.png b/tutorial/figures/fig4.png index 09c6593..f047f60 100644 Binary files a/tutorial/figures/fig4.png and b/tutorial/figures/fig4.png differ diff --git a/tutorial/figures/fig_model.png b/tutorial/figures/fig_model.png index 933afd6..6d0f365 100644 Binary files a/tutorial/figures/fig_model.png and b/tutorial/figures/fig_model.png differ diff --git a/tutorial/tutorial.jl b/tutorial/tutorial.jl index 0530a32..49cdda2 100644 --- a/tutorial/tutorial.jl +++ b/tutorial/tutorial.jl @@ -15,19 +15,25 @@ #md # these neurotransmitters change over time, and what the new steady-state concentrations are. # ## 2. Key concepts of Glutamate and GABA cycling -#md # Glutamate and GABA are the primary excitatory and inhibitory neurotransmitters in the brain, -#md # synthesized de novo in neurons, with GABA synthesis being specific to GABAergic neurons. -#md # Due to the high energetic cost of neurotransmission, approximately 80% of released glutamate -#md # and GABA is recycled by astrocytes via glutamine. This model incorporates the essential steps of -#md # glutamate and GABA cycling, using Michaelis–Menten kinetics, to simulate the cycling rate -#md # and the resultant concentrations. - -#md # Another key aspect of neurotransmitter cycling is its coupling to glucose oxidation, -#md # a phenomenon widely observed in experiments. -#md # Recent studies suggest that this coupling links glucose oxidation to the conversion -#md # of glutamine to glutamate in neurons. -#md # This process has been described as the pseudo–malate–aspartate shuttle (PMAS). -#md # In this model, PMAS exerts strong flux control over the cycling rate. +#md # Glutamate and GABA are the primary excitatory and inhibitory neurotransmitters in the brain. +#md # Both are synthesized de novo in neurons, with GABA synthesis being specific to +#md # GABAergic neurons. +#md # +#md # Following synaptic release during neurotransmission, approximately 80% of released +#md # glutamate and GABA is taken up by astrocytes and converted to glutamine, which is +#md # subsequently transported back to neurons and converted to glutamate and GABA. +#md # The recycling process is closely coupled to neuronal energy metabolism, helping +#md # maintain a balance between neuronal signaling and its energetic demands across +#md # varying levels of activity. One proposed coupling mechanism is the +#md # pseudo-malate-aspartate shuttle (PMAS). +#md # +#md # The current model incorporates the key enzymatic reactions involved in +#md # neurotransmitter cycling in neurons and astrocytes and couples them to neuronal +#md # glucose metabolism via the PMAS. It integrates algebraic expressions and +#md # differential equations into a unified framework, with parameters derived primarily +#md # from in vivo studies. In addition, the model includes equations describing ketone +#md # pharmacokinetics in the brain, enabling the simulation of ketosis and its effects +#md # on neurotransmitter cycling. #md # ![](figures/fig_model.png) @@ -55,7 +61,7 @@ prob = ODEProblem(ntmodel.model, ntmodel.u0, (0.0, 100.0), ntmodel.ps); sol = solve(prob, Tsit5()); # Then plot the results: -fig = Figure(resolution=(400, 300)) +fig = Figure(size=(400, 300)) ax = Axis(fig[1, 1]) for state in [:GLU_e, :GABA_i, :GLN_a] lines!(ax, sol.t, sol[state], linewidth=2, label=string(state)) @@ -80,7 +86,7 @@ prob = ODEProblem(ntmodel.model, ntmodel.u0, (0.0, 300.0), ntmodel.ps); sol = solve(prob, Tsit5()); #md # And then, we'll check how the system responds to the change in initial conditions: -fig = Figure(resolution=(400, 300)) +fig = Figure(size=(400, 300)) ax = Axis(fig[1, 1]) for state in [:GLU_e, :GABA_i, :GLN_a] lines!(ax, sol.t, sol[state], linewidth=2, label=string(state)) @@ -114,7 +120,7 @@ prob = ODEProblem(ntmodel.model, ntmodel.u0, (0.0, 300.0), ntmodel.ps); sol = solve(prob, Tsit5()); # Plot the results: -fig = Figure(resolution=(400, 300)) +fig = Figure(size=(400, 300)) ax = Axis(fig[1, 1]) for state in [:GLU_e, :GABA_i, :GLN_a] lines!(ax, sol.t, sol[state], linewidth=2, label=string(state)) @@ -129,11 +135,18 @@ display(fig) #md # We can see that in this case, the system returns to a new steady state after the perturbation. # ### Example 2: Effect of GABA-T inhibition -#md # In this example, we will simulate the effect of GABA-T inhibition. -#md # GABA-T is an enzyme that catalyzes the conversion of GABA to Glutamate in GABAergic neurons. -#md # We will simulate the effects of pharmacological inhibition of GABA-T by reducing the Vmax parameter -#md # of this enzyme. The corresponding parameter is labeled as Vmax_GT in the model. In this scenario, we -#md # will simulate for longer due to slower dynamics. +#md # In this example, we will simulate the effect of ketosis on neurotransmitter cycling. +#md # Ketosis is a metabolic state characterized by elevated levels of ketone bodies in the +#md # blood. +#md # Ketone bodies are transported from the blood into the brain, where they can serve as +#md # an alternative energy substrate alongside glucose, thereby altering neuronal +#md # metabolism. This may consequently affect neurotransmitter cycling, as it is tightly +#md # coupled to cellular energy metabolism. +#md # Here, we use the model to predict how ketosis affects neuronal metabolism and +#md # neurotransmitter cycling. +#md # The model parameter corresponding to plasma ketone concentration is KetP. +#md # To simulate ketosis, we will set the initial ketone concentration to 3 mM +#md # (originally 0 mM) and run the simulation. # Start with a fresh instance of the NTModel: ntmodel = NTModel(); @@ -147,7 +160,7 @@ sol = solve(prob, Tsit5()); #md # We will plot GABA_i only as the Glutamate/Glutamine concentrations are #md # not expected to change significantly: -fig = Figure(resolution=(400, 300)) +fig = Figure(size=(400, 300)) ax = Axis(fig[1, 1]) lines!(ax, sol.t, sol[:GABA_i], linewidth=2, label="GABA_i") axislegend(ax, position=:rb) diff --git a/tutorial/tutorial.md b/tutorial/tutorial.md index e4460c1..2aee52c 100644 --- a/tutorial/tutorial.md +++ b/tutorial/tutorial.md @@ -1,3 +1,7 @@ +```@meta +EditURL = "tutorial.jl" +``` + # Neuroblox Tutorial: Neurotransmitter Cycling Model This tutorial will show you how to use the Neurotransmitter cycling module of Neuroblox. @@ -14,19 +18,25 @@ or pharmacological interventions. The simulations will output how the concentrat these neurotransmitters change over time, and what the new steady-state concentrations are. ## 2. Key concepts of Glutamate and GABA cycling -Glutamate and GABA are the primary excitatory and inhibitory neurotransmitters in the brain, -synthesized de novo in neurons, with GABA synthesis being specific to GABAergic neurons. -Due to the high energetic cost of neurotransmission, approximately 80% of released glutamate -and GABA is recycled by astrocytes via glutamine. This model incorporates the essential steps of -glutamate and GABA cycling, using Michaelis–Menten kinetics, to simulate the cycling rate -and the resultant concentrations. - -Another key aspect of neurotransmitter cycling is its coupling to glucose oxidation, -a phenomenon widely observed in experiments. -Recent studies suggest that this coupling links glucose oxidation to the conversion -of glutamine to glutamate in neurons. -This process has been described as the pseudo–malate–aspartate shuttle (PMAS). -In this model, PMAS exerts strong flux control over the cycling rate. +Glutamate and GABA are the primary excitatory and inhibitory neurotransmitters in the brain. +Both are synthesized de novo in neurons, with GABA synthesis being specific to +GABAergic neurons. + +Following synaptic release during neurotransmission, approximately 80% of released +glutamate and GABA is taken up by astrocytes and converted to glutamine, which is +subsequently transported back to neurons and converted to glutamate and GABA. +The recycling process is closely coupled to neuronal energy metabolism, helping +maintain a balance between neuronal signaling and its energetic demands across +varying levels of activity. One proposed coupling mechanism is the +pseudo-malate-aspartate shuttle (PMAS). + +The current model incorporates the key enzymatic reactions involved in +neurotransmitter cycling in neurons and astrocytes and couples them to neuronal +glucose metabolism via the PMAS. It integrates algebraic expressions and +differential equations into a unified framework, with parameters derived primarily +from in vivo studies. In addition, the model includes equations describing ketone +pharmacokinetics in the brain, enabling the simulation of ketosis and its effects +on neurotransmitter cycling. ![](figures/fig_model.png) @@ -37,7 +47,7 @@ First, make sure you have the module file in your workspace, then follow the ste Import the necessary packages: -````julia +````@example tutorial using OrdinaryDiffEq using Catalyst using CairoMakie @@ -46,8 +56,9 @@ using Main.NTModelModule Create an instance of the NTModel: -````julia +````@example tutorial ntmodel = NTModel(); +nothing #hide ```` By default, the model is set up with initial conditions that correspond to a steady state. @@ -56,15 +67,16 @@ The parameters are defined in the `ps` field. From the model, we can create an ODEProblem as usual and solve it: -````julia +````@example tutorial prob = ODEProblem(ntmodel.model, ntmodel.u0, (0.0, 100.0), ntmodel.ps); sol = solve(prob, Tsit5()); +nothing #hide ```` Then plot the results: -````julia -fig = Figure(resolution=(400, 300)) +````@example tutorial +fig = Figure(size=(400, 300)) ax = Axis(fig[1, 1]) for state in [:GLU_e, :GABA_i, :GLN_a] lines!(ax, sol.t, sol[state], linewidth=2, label=string(state)) @@ -83,22 +95,24 @@ Now let's change one of the initial conditions to see how the model responds. We'll save the original value of glutamate in excitatory neurons (10 mM) and then set it to 12 mM: -````julia +````@example tutorial GLU_e_original = ntmodel.u0[:GLU_e]; ntmodel.u0[:GLU_e] = 12.0; +nothing #hide ```` We'll run the simulation again with the new initial conditions: -````julia +````@example tutorial prob = ODEProblem(ntmodel.model, ntmodel.u0, (0.0, 300.0), ntmodel.ps); sol = solve(prob, Tsit5()); +nothing #hide ```` And then, we'll check how the system responds to the change in initial conditions: -````julia -fig = Figure(resolution=(400, 300)) +````@example tutorial +fig = Figure(size=(400, 300)) ax = Axis(fig[1, 1]) for state in [:GLU_e, :GABA_i, :GLN_a] lines!(ax, sol.t, sol[state], linewidth=2, label=string(state)) @@ -124,27 +138,30 @@ and run the simulation. Start with a fresh instance of the NTModel: -````julia +````@example tutorial ntmodel = NTModel(); +nothing #hide ```` Set the initial concentration of plasma ketone bodies to 3 mM: -````julia +````@example tutorial ntmodel.ps[:KET_p] = 3.0; +nothing #hide ```` Run the simulation: -````julia +````@example tutorial prob = ODEProblem(ntmodel.model, ntmodel.u0, (0.0, 300.0), ntmodel.ps); sol = solve(prob, Tsit5()); +nothing #hide ```` Plot the results: -````julia -fig = Figure(resolution=(400, 300)) +````@example tutorial +fig = Figure(size=(400, 300)) ax = Axis(fig[1, 1]) for state in [:GLU_e, :GABA_i, :GLN_a] lines!(ax, sol.t, sol[state], linewidth=2, label=string(state)) @@ -160,36 +177,46 @@ display(fig) We can see that in this case, the system returns to a new steady state after the perturbation. ### Example 2: Effect of GABA-T inhibition -In this example, we will simulate the effect of GABA-T inhibition. -GABA-T is an enzyme that catalyzes the conversion of GABA to Glutamate in GABAergic neurons. -We will simulate the effects of pharmacological inhibition of GABA-T by reducing the Vmax parameter -of this enzyme. The corresponding parameter is labeled as Vmax_GT in the model. In this scenario, we -will simulate for longer due to slower dynamics. +In this example, we will simulate the effect of ketosis on neurotransmitter cycling. +Ketosis is a metabolic state characterized by elevated levels of ketone bodies in the +blood. +Ketone bodies are transported from the blood into the brain, where they can serve as +an alternative energy substrate alongside glucose, thereby altering neuronal +metabolism. This may consequently affect neurotransmitter cycling, as it is tightly +coupled to cellular energy metabolism. +Here, we use the model to predict how ketosis affects neuronal metabolism and +neurotransmitter cycling. +The model parameter corresponding to plasma ketone concentration is KetP. +To simulate ketosis, we will set the initial ketone concentration to 3 mM +(originally 0 mM) and run the simulation. Start with a fresh instance of the NTModel: -````julia +````@example tutorial ntmodel = NTModel(); +nothing #hide ```` Reduce the Vmax parameter of GABA-T by 5%: -````julia +````@example tutorial ntmodel.ps[:Vmax_GT] *= 0.95; +nothing #hide ```` Run the simulation -````julia +````@example tutorial prob = ODEProblem(ntmodel.model, ntmodel.u0, (0.0, 7200.0), ntmodel.ps); sol = solve(prob, Tsit5()); +nothing #hide ```` We will plot GABA_i only as the Glutamate/Glutamine concentrations are not expected to change significantly: -````julia -fig = Figure(resolution=(400, 300)) +````@example tutorial +fig = Figure(size=(400, 300)) ax = Axis(fig[1, 1]) lines!(ax, sol.t, sol[:GABA_i], linewidth=2, label="GABA_i") axislegend(ax, position=:rb)