# Mitochondrial Dynamics

Thesis: Mathematical Modeling and Computational Simulation of

Metabolic Regulation of Mitochondrial Dynamics in Aging by Chen Chang

## Fitting AMPK activity

The numbers come from Peter's thesis.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
# model
def ampk_activity(r, vmax, k_ampk, n):
    keq_adk = 0.5
    atp = 1
    adp = r
    amp = adp**2 * keq_adk / atp
    # Normalize
    amp = amp / (amp + adp + atp)
    return vmax * (amp**n) / (amp**n + k_ampk**n)

# bounds
lb = [0.0, 1e-5]
ub = [35.0, 1.0]
# data
xdata = np.array([0.03355, 0.037518, 0.04615, 0.081028, 0.203052])
ydata = np.array([10.0, 15.0, 20.0, 25.0, 30.0])
def ampk_activity_n1(r, vmax, k_ampk):
    return ampk_activity(r, vmax, k_ampk, 1.0)
popt, pcov = curve_fit(ampk_activity_n1, xdata, ydata, bounds=(lb, ub))
print('vmax, k_ampk =', popt)
vmax, k_ampk = popt
xs = np.linspace(0, 0.2)
ys = ampk_activity(xs, vmax, k_ampk, 1.0)
fig, ax = plt.subplots()
ax.plot(xs, ys)
ax.scatter(xdata, ydata, c='r')
fig
Python

### Considering cooperativity

# bounds
lb = [30.0, 0.0, 1.0]
ub = [35.0, 1.0, 5.0]
popt, pcov = curve_fit(ampk_activity, xdata, ydata, bounds=(lb, ub))
print('vmax, k_ampk, n =', popt)
vmax, k_ampk, n = popt
xs = np.linspace(0, 0.2)
ys = ampk_activity(xs, vmax, k_ampk, n)
fig, ax = plt.subplots()
ax.plot(xs, ys)
ax.scatter(xdata, ydata, c='r')
fig
Python

## The mitochondrion model

Enhanced BLPS model, with the addition of adenylate kinase (AK) to perform conversions between ATP, ADP, and AMP. The model behavior is highly sensitive to MCUand ANT maximal rates.

# enhanced BLPS (eBLPS) model for mitochondrial ATP production in a minimistic way
using Parameters, LabelledArrays
"Converting energy charge (EC) to fractions of adenylates (ATP, ADP, AMP)"
function energy_charge_to_adenylates(ec, keq = 1)
    k = 8 * keq - 2
    t = ec + (1 - sqrt(1 + 2 * k * ec * (1 - ec))) / k
    d = 2 * (ec - t)
    m = 1 - t - d
    return (t, d, m)
end
# Faraday constant (sA/mol)
const F = 96485
# Ideal gas constant
const R = 8.314
# Temperature (K)
const T = 310
# Inerse of specific Nernst voltage
const F_RT = F / (R * T)
# Baseline values
const NAD0 = 10E-3   # M
const A0 = 15E-3	 # M
const CA0 = 0.2E-6	 # M
const ΔΨ0 = 0.091    # V
const FBP0 = 1E-6    # M
# Boltzman factor for a specific voltage
pv(v) = exp(-x * F_RT)
# Michaelis-Menten
_mm(x, k = one(x)) = x / (x + k)
# Hill function
_hill(x, k, n) = _mm(x^n, k^n)
# Sigmoid function
σ(x, a = one(x)) = _mm(one(x), a * exp(-x))
"Pyruvate dehydrogenase (PDH) parameters"
@with_kw struct PDHParams
    V0 = 0.2E-3         # M/s
    FBP = 2E-6          # M
    VMAX = V0 * sqrt(FBP / FBP0)  # M/s
    K_CA = 0.05 * CA0   # M
    K_NAD = 1		    # dimless
end
"PDH rate"
function jpdh(nadh, nad, ca_m, p::PDHParams)
    @unpack VMAX, K_CA, K_NAD = p
    jPDH = VMAX * _mm(ca_m, K_CA) * nad / (K_NAD * nad + nadh)
end
"Electron transport chain parameters"
@with_kw struct ETCParams
    V0 = 0.6E-3              # M/s
    K_NADH = 0.01 * NAD0     # M
    A4 = 4.23E-16			 # dimless
    A5 = 18.2 / ΔΨ0          # 1/V
    A14 = 11.7				 # dimless
end
"Oxygen consumption rate"
function jo(nadh, ΔΨ, p::ETCParams)
    @unpack V0, K_NADH, A4, A5 = p
    jO = V0 * _mm(nadh, K_NADH) * σ(-A5 * ΔΨ, A4)
end
"Proton flux by ETC"
jhres(jO, p::ETCParams) = jO * p.A14
"ATP synthase parameters"
@with_kw struct F1Params
    V0 = 23.3E-3 / 0.67      # M/s
    A12 = 5.1E9              # dimless
    A13 = 10.7 / ΔΨ0		 # 1/V
    K_ATP = 0.67 * A0        # M
    A15 = 3.43			     # dimless
end
"ATP synthase rate for ATP"
function jf1fo(atp_m, ΔΨ, p::F1Params)
    @unpack V0, K_ATP, A12, A13 = p
    jF1Fo = V0 * _mm(K_ATP, atp_m) * σ(A13 * ΔΨ, A12)
end
"Proton flux through ATP synthase"
jhatp(jF1Fo, p::F1Params) = jF1Fo * p.A15
"Protn leak parameters"
@with_kw struct HLeakParams
    V0 = 0.182e-3 / ΔΨ0      # M/Vs
    Δ = 0.16 * ΔΨ0			 # V
end
"Proton leak rate (linear)"
jhleak(ΔΨ, p::HLeakParams) = p.V0 * (ΔΨ - p.Δ)
"Adenylate Nucleotide Translocase (ANT) parameters"
@with_kw struct ANTParams
    V0 = 5E-3 * 0.2      # M/s
    A6 = 1 - 0.861       # dimless
    A7 = 0.139			 # dimless
    A8 = 0.111			 # dimless
    A9 = 3.37 / ΔΨ0	     # 1/V
    A10 = 1.68 / ΔΨ0     # 1/V
end
"ANT rates"
function jant(atp_m, adp_m, ΔΨ, p::ANTParams)
    @unpack V0, A6, A7, A8, A9, A10 = p
    a = atp_m / (adp_m + A6 * atp_m)
    b = adp_m * (A7 - A8 * exp(-A9 * ΔΨ))
    c = adp_m + A8 * atp_m * exp(-A10 * ΔΨ)
    jANT = V0 * a * b / c
end
"mitochondrial NCX (NCLX) parameters"
@with_kw struct NCLXParams
    V0 = 0.1E-5				    # M/s
    A17 = 1.46 / ΔΨ0			# 1/V
end
"NCLX rate"
jnaca(ca_m, ca_i, ΔΨ, p::NCLXParams) = p.V0 * ca_m / ca_i * exp(p.A17 * ΔΨ)
"mitochondrial Calcium uniporter (MCU) parameters"
@with_kw struct MCUParams
    V0 = 0.11e-3 / 0.01 * 0.1 # M/s
    A18 = 6.73 / ΔΨ0        # 1/V
    A19 = 0.52 / CA0	    # 1/M
    A20 = 0.01 / CA0        # 1/M
    L = 110  				# dimless
    NA = 2.8				# dimless
end
"MCU rate"
function juni(ca_i, ΔΨ, p::MCUParams)
    @unpack A18 = p
    δ = -A18 * (ΔΨ - ΔΨ0)
    fV = δ / expm1(δ)
    @unpack A19, A20, L, NA, V0 = p
    relax = A20 * ca_i
    taut = (1 + A19 * ca_i)^NA * (1 + relax)^3
    fMWC = relax * taut / (L + taut * (1 + relax))
    jUni = V0 * fV * fMWC
end
"Calcium wave parameters"
@with_kw struct CaWaveParams
    BASELINE = 1.2E-6
    A = 0.1E-6
    PERIOD = 150
    FREQ_2 = 2 * inv(PERIOD)
end
"Intracellular calcium level"
cai(t, p::CaWaveParams) = p.BASELINE + p.A * sinpi(p.FREQ_2 * t)
"GTPase parameters"
@with_kw struct GGCLPParams
    VMAX = 0.58E-3   # M/s (Modified by 100x)
    KGTP = 0.12E-3	 # M (Modified by 100x)
    ΣGTP = 1.5E-3    # Total GTP pool (M)
    R_ANT = 0.1      # Flux of ANT goes to GTP (dimless)
end
"GTPase rate"
jggclp(gtp, p::GGCLPParams) = p.VMAX * _mm(gtp, p.KGTP)
"GTP production rates"
jgtp(gtp, jANT, p::GGCLPParams) = p.R_ANT * jANT * (1 - gtp / p.ΣGTP)
"Adenylate kinase (ADK) parameters"
@with_kw struct ADKParams
    KF = 0.73         # M/s
    KEQ = 1.0		  # dimless
    KB = KF / KEQ
end
"ADK rates"
jadk(atp, adp, amp, p::ADKParams) = p.KF * adp^2 - p.KB * (atp * amp)
"Get ADP concentration by applying Quasi-steady-state assumption (QSSA) on AMP"
function adp_qssa(atp, ΣA, p::ADKParams)
    ke = p.KEQ
    adp = (sqrt(atp^2 + 4 * ke * atp * (ΣA - atp)) - atp) / (2 * ke)
end
"Get AMP concentration by applying Quasi-steady-state assumption (QSSA) on AMP"
function amp_qssa(atp, ΣA, p::ADKParams)
    adp = adp_qssa(atp, p)
    amp = ΣA - atp - adp
end
@with_kw struct BLPSParams
    # NAD+ + NADH concentration (M)
    ΣNAD = 10e-3
    # ATP + ADP + AMP concentration in the mitochodria (M)
    ΣA = 10e-3
    # IMM capacitance (M/V)
    CM = 1.8e-3
    CM_i = inv(CM)
    # Calcium buffering factor
    δCa = 0.01
    # Parameters from other components
    pPDH = PDHParams()
    pETC = ETCParams()
    pF1 = F1Params()
    pLeak = HLeakParams()
    pANT = ANTParams()
    pNCLX = NCLXParams()
    pMCU = MCUParams()
    pCai = CaWaveParams()
    pGGCLP = GGCLPParams()
    pADK = ADKParams()
end
"Calculate initial conditions"
function get_u0(p::BLPSParams; nadhRatio = 0.01, energyCharge = 0.75, ca_m = 0.6e-6, ΔΨ = 0.0, gtp = 0.0)
    @unpack ΣNAD, ΣA = p
    t, _, _ = energy_charge_to_adenylates(energyCharge, p.pADK.KEQ)
    u0 = LVector(
        nadh = ΣNAD * nadhRatio,
        atp_m = t * ΣA,
        ca_m = ca_m,
        ΔΨ = ΔΨ,
        gtp = gtp
        )
end
"enhanced BLPS model"
function blps!(du, u, p::BLPSParams, t)
    @unpack ΣNAD, ΣA = p
    nadh, atp_m, ca_m, ΔΨ, gtp = u.nadh, u.atp_m, u.ca_m, u.ΔΨ, u.gtp
    nad = ΣNAD - nadh
    @unpack pPDH, pCai, pETC, pF1, pLeak, pANT, pNCLX, pMCU, pGGCLP, pADK = p
    adp_m = adp_qssa(atp_m, ΣA, pADK)
    ca_i = cai(t, pCai)
    jPDH = jpdh(nadh, nad, ca_m, pPDH)
    jO = jo(nadh, ΔΨ, pETC)
    jHRes = jhres(jO, pETC)
    jF1 = jf1fo(atp_m, ΔΨ, pF1)
    jHATP = jhatp(jF1, pF1)
    jHleak = jhleak(ΔΨ, pLeak)
    jANT = jant(atp_m, adp_m, ΔΨ, pANT)
    jNaCa = jnaca(ca_m, ca_i, ΔΨ, pNCLX)
    jUNi = juni(ca_i, ΔΨ, pMCU)
    jGTP = jgtp(gtp, jANT, pGGCLP)
    jGGCLP =  jggclp(gtp, pGGCLP)
    @unpack CM_i, δCa = p
    du.nadh = jPDH - jO
    du.atp_m = jF1 - jANT
    du.ca_m = δCa * (jUNi - jNaCa)
    du.ΔΨ = CM_i * (jHRes - jANT - jNaCa - 2 * jUNi - jHATP - jHleak)
    du.gtp = jGTP - jGGCLP
    return (jPDH = jPDH,
            jO = jO,
            jHRes = jHRes,
            jF1 = jF1,
            jHATP = jHATP,
            jHleak = jHleak,
            jANT = jANT,
            jNaCa = jNaCa,
            jUNi = jUNi,
            jGTP = jGTP,
            jGGCLP = jGGCLP)
end
4.4s
blps!
using DifferentialEquations, Plots, BenchmarkTools
# Plots energy charge relationship
p = plot(title = "Energy charge (Keq = 1.0)", legend = :top, xlabel = "EC", ylabel = "Fraction")
for (i, la) in enumerate(["ATP", "ADP", "AMP"])
    plot!(p, x -> energy_charge_to_adenylates(x, 1.0)[i], 0.0, 1.0, label = la)
end
display(p)
74.5s
# Parameters for quasi-steady state solution
# Tuned down ANT and MCU activities
p = BLPSParams()
u0 = get_u0(p; nadhRatio = 0.0088, energyCharge = 0.9, ca_m = 0.6e-6, ΔΨ = 0.14, gtp = 0.0)
tspan = (0.0, 10000.0)
prob = ODEProblem(blps!, u0, tspan, p)
sol = solve(prob)
65.7s
plot(sol, vars = (0, 1), label = "NADH", lw = 2)
7.1s
plot(sol, vars = (0, 2), label = "ATP", legend=:right)
2.3s
plot(sol, vars = (0, 3), label = "Ca(mito)", legend = :right, lw=1)
plot!(t->cai(t, p.pCai), label = "Ca(cyto)", tspan[1], tspan[2], legend = :right)
5.0s
plot(sol, vars = (0, 4), label = "psi", lw = 1)
2.2s
plot(sol, vars = (0, 5), label = "GTP", lw = 2)
2.2s
# Aged mitochondria
# Increased frequency of Ca wave (100sec vs. 150sec)
# Increased ATP synthase activity (1.47x)
# Increased proton leak (4x)
pAged = BLPSParams(
   p;
   pCai = CaWaveParams(PERIOD = 100),
   pF1 = F1Params(V0 = p.pF1.V0 * 1.47),
   pLeak = HLeakParams(V0 = p.pLeak.V0 * 4)
)
u0Aged = get_u0(pAged; nadhRatio = 0.01, energyCharge = 0.75, ca_m = 0.6e-6, ΔΨ = 0.0, gtp = 0.0)
tspan = (0.0, 1000.0)
probAged = ODEProblem(blps!, u0Aged, tspan, pAged)
solAged = solve(probAged, dt = 1e-3)
2.1s
plot(solAged, vars = (0, 1), label = "NADH (Aged)", ylim = (8.75e-5, 9.0e-5))
1.4s
plot(solAged, vars = (0, 2), label = "ATP (Aged)")
1.0s
plot(solAged, vars = (0, 3), label = "Cam(Aged)")
0.8s
plot(solAged, vars = (0, 4), label = "psi(Aged)", ylims = (0.14, 0.16), lw = 1)
0.9s
plot(solAged, vars = (0, 5), label = "GTP(Aged)")
1.0s

## The AMPK activation model

https://nextjournal.com/sosiristseng/fitting-ampk-activity/

0.2s