Metabolism Simulation Programs

Interactive simulations exploring the quantitative and computational side of human metabolism. From enzyme kinetics and glycolytic flux to whole-body glucose homeostasis and cancer metabolism, these programs let you adjust parameters and observe how metabolic systems respond in real time.

1. Michaelis-Menten Enzyme Kinetics

The Michaelis-Menten equation is the cornerstone of enzyme kinetics, describing how reaction velocity depends on substrate concentration. For a simple enzyme-catalyzed reaction where enzyme E binds substrate S to form product P via an intermediate complex ES, the steady-state rate equation is:

$$v = \frac{V_{\max}[S]}{K_m + [S]}$$

Here $V_{\max} = k_{\text{cat}}[E]_T$ is the maximum velocity when all enzyme is saturated, and $K_m$ is the Michaelis constant, equal to the substrate concentration at which $v = V_{\max}/2$. A low $K_m$ indicates high substrate affinity. The Lineweaver-Burk double-reciprocal plot linearizes this relationship:

$$\frac{1}{v} = \frac{K_m}{V_{\max}} \cdot \frac{1}{[S]} + \frac{1}{V_{\max}}$$

This simulation plots both the hyperbolic Michaelis-Menten curve and its Lineweaver-Burk linearization, annotating the key kinetic parameters graphically. Adjust the parameters below to explore how enzyme affinity and catalytic capacity shape the velocity profile.

Python

Parameters

Vmax (µmol/min)

Maximum velocity

Km (mM)

Michaelis constant

S range max (mM)

Maximum substrate concentration

script.py90 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Parameters
VMAX_VALUE = 100.0
KM_VALUE = 5.0
SMAX_VALUE = 50.0

Vmax = VMAX_VALUE
Km = KM_VALUE
S_max = SMAX_VALUE

S = np.linspace(0.001, S_max, 500)
v = Vmax * S / (Km + S)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
fig.patch.set_facecolor('#0f172a')

# --- Left: Michaelis-Menten curve ---
ax1.set_facecolor('#1e293b')
ax1.plot(S, v, color='#4ade80', linewidth=2.5, label='v = Vmax[S]/(Km+[S])')
ax1.axhline(y=Vmax, color='#f87171', linestyle='--', linewidth=1.5, label=f'Vmax = {Vmax:.1f}')
ax1.axhline(y=Vmax/2, color='#fbbf24', linestyle=':', linewidth=1.2, alpha=0.7, label=f'Vmax/2 = {Vmax/2:.1f}')
ax1.axvline(x=Km, color='#60a5fa', linestyle=':', linewidth=1.2, alpha=0.7, label=f'Km = {Km:.1f} mM')
ax1.plot(Km, Vmax/2, 'o', color='#fbbf24', markersize=10, zorder=5)
ax1.annotate(f'(Km, Vmax/2)\n({Km:.1f}, {Vmax/2:.1f})',
             xy=(Km, Vmax/2), xytext=(Km + S_max*0.1, Vmax*0.35),
             color='#fbbf24', fontsize=10,
             arrowprops=dict(arrowstyle='->', color='#fbbf24', lw=1.5))
ax1.set_xlabel('[S] (mM)', fontsize=12)
ax1.set_ylabel('v (umol/min)', fontsize=12)
ax1.set_title('Michaelis-Menten Kinetics', fontsize=14)
ax1.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1', fontsize=9)
ax1.tick_params(colors='#cbd5e1')
ax1.xaxis.label.set_color('#cbd5e1')
ax1.yaxis.label.set_color('#cbd5e1')
ax1.title.set_color('#e2e8f0')
for spine in ax1.spines.values():
    spine.set_color('#334155')
ax1.set_xlim(0, S_max)
ax1.set_ylim(0, Vmax * 1.15)
ax1.grid(True, alpha=0.15, color='#cbd5e1')

# --- Right: Lineweaver-Burk ---
ax2.set_facecolor('#1e293b')
S_lb = np.linspace(Km * 0.15, S_max, 300)
v_lb = Vmax * S_lb / (Km + S_lb)
inv_S = 1.0 / S_lb
inv_v = 1.0 / v_lb

ax2.plot(inv_S, inv_v, color='#c084fc', linewidth=2.5, label='1/v vs 1/[S]')

# Extend line to y-axis intercept
slope = Km / Vmax
y_int = 1.0 / Vmax
x_int = -1.0 / Km
x_range = np.linspace(x_int * 0.8, max(inv_S) * 1.05, 100)
y_range = slope * x_range + y_int
ax2.plot(x_range, y_range, '--', color='#c084fc', alpha=0.4, linewidth=1.5)

ax2.axhline(y=0, color='#475569', linewidth=0.8)
ax2.axvline(x=0, color='#475569', linewidth=0.8)
ax2.plot(0, y_int, 's', color='#f87171', markersize=10, zorder=5, label=f'y-int = 1/Vmax = {y_int:.4f}')
ax2.plot(x_int, 0, 'D', color='#60a5fa', markersize=10, zorder=5, label=f'x-int = -1/Km = {x_int:.3f}')

ax2.set_xlabel('1/[S] (1/mM)', fontsize=12)
ax2.set_ylabel('1/v (min/umol)', fontsize=12)
ax2.set_title('Lineweaver-Burk Plot', fontsize=14)
ax2.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1', fontsize=9)
ax2.tick_params(colors='#cbd5e1')
ax2.xaxis.label.set_color('#cbd5e1')
ax2.yaxis.label.set_color('#cbd5e1')
ax2.title.set_color('#e2e8f0')
for spine in ax2.spines.values():
    spine.set_color('#334155')
ax2.grid(True, alpha=0.15, color='#cbd5e1')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a', edgecolor='none')

print(f"Michaelis-Menten Kinetics Analysis")
print(f"={'='*40}")
print(f"Vmax = {Vmax:.2f} umol/min")
print(f"Km   = {Km:.2f} mM")
print(f"Catalytic efficiency (Vmax/Km) = {Vmax/Km:.2f} umol/(min*mM)")
print(f"At [S] = Km: v = {Vmax*Km/(Km+Km):.2f} umol/min (= Vmax/2)")
print(f"At [S] = 10*Km: v = {Vmax*10*Km/(Km+10*Km):.2f} umol/min ({Vmax*10*Km/(Km+10*Km)/Vmax*100:.1f}% of Vmax)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Understanding the Results

The left plot shows the classic hyperbolic saturation curve. At low $[S] \ll K_m$, the rate increases approximately linearly with slope $V_{\max}/K_m$ (catalytic efficiency). At high $[S] \gg K_m$, the enzyme is saturated and the rate plateaus at $V_{\max}$. The yellow dot marks the half-saturation point.

The Lineweaver-Burk plot is the double-reciprocal linearization. The y-intercept gives $1/V_{\max}$ and the x-intercept gives $-1/K_m$. While useful historically for parameter estimation, this transform amplifies errors at low substrate concentrations.

2. Enzyme Inhibition Comparison

Enzyme inhibitors are central to pharmacology and metabolic regulation. The four main types of reversible inhibition each produce characteristic changes in the apparent kinetic parameters. Competitive inhibitors bind the active site, increasing the apparent $K_m$ without affecting $V_{\max}$:

$$v_{\text{competitive}} = \frac{V_{\max}[S]}{K_m\left(1 + \frac{[I]}{K_i}\right) + [S]}$$

Uncompetitive inhibitors bind only the ES complex, decreasing both apparent $K_m$ and apparent $V_{\max}$ by the same factor $\alpha' = 1 + [I]/K_i'$. Mixed (noncompetitive when $K_i = K_i'$) inhibitors bind both E and ES, reducing$V_{\max}$ while changing $K_m$ depending on relative binding affinities.

$$v_{\text{mixed}} = \frac{V_{\max}[S]}{\alpha K_m + \alpha' [S]}, \quad \alpha = 1+\frac{[I]}{K_i}, \quad \alpha' = 1+\frac{[I]}{K_i'}$$

Python

Parameters

Vmax (µmol/min)

Maximum velocity

Km (mM)

Michaelis constant

Ki (mM)

Inhibition constant

[I] (mM)

Inhibitor concentration

script.py96 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

VMAX_VALUE = 100.0
KM_VALUE = 5.0
KI_VALUE = 3.0
INHIB_VALUE = 5.0

Vmax = VMAX_VALUE
Km = KM_VALUE
Ki = KI_VALUE
I_conc = INHIB_VALUE

S = np.linspace(0.01, 60, 500)

# No inhibitor
v_none = Vmax * S / (Km + S)

# Competitive: apparent Km increases
alpha_c = 1 + I_conc / Ki
v_comp = Vmax * S / (Km * alpha_c + S)

# Uncompetitive: both apparent Km and Vmax decrease
alpha_u = 1 + I_conc / Ki
v_uncomp = (Vmax / alpha_u) * S / (Km / alpha_u + S)

# Mixed (noncompetitive, Ki = Ki')
alpha_m = 1 + I_conc / Ki
v_mixed = (Vmax / alpha_m) * S / (Km + S)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
fig.patch.set_facecolor('#0f172a')

colors = ['#4ade80', '#f87171', '#60a5fa', '#fbbf24']
labels = ['No inhibitor', 'Competitive', 'Uncompetitive', 'Mixed (noncompetitive)']
curves = [v_none, v_comp, v_uncomp, v_mixed]

# Left: v vs [S]
ax1.set_facecolor('#1e293b')
for c, lbl, clr in zip(curves, labels, colors):
    ax1.plot(S, c, linewidth=2.2, color=clr, label=lbl)
ax1.set_xlabel('[S] (mM)', fontsize=12)
ax1.set_ylabel('v (umol/min)', fontsize=12)
ax1.set_title('Inhibition: v vs [S]', fontsize=14)
ax1.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1', fontsize=9)
ax1.tick_params(colors='#cbd5e1')
ax1.xaxis.label.set_color('#cbd5e1')
ax1.yaxis.label.set_color('#cbd5e1')
ax1.title.set_color('#e2e8f0')
for spine in ax1.spines.values():
    spine.set_color('#334155')
ax1.grid(True, alpha=0.15, color='#cbd5e1')
ax1.set_xlim(0, 60)
ax1.set_ylim(0, Vmax * 1.15)

# Right: Lineweaver-Burk
ax2.set_facecolor('#1e293b')
S_lb = np.linspace(1.0, 60, 300)
for c_func, lbl, clr in zip(
    [lambda s: Vmax * s / (Km + s),
     lambda s: Vmax * s / (Km * alpha_c + s),
     lambda s: (Vmax / alpha_u) * s / (Km / alpha_u + s),
     lambda s: (Vmax / alpha_m) * s / (Km + s)],
    labels, colors):
    v_lb = c_func(S_lb)
    ax2.plot(1/S_lb, 1/v_lb, linewidth=2.2, color=clr, label=lbl)

ax2.set_xlabel('1/[S] (1/mM)', fontsize=12)
ax2.set_ylabel('1/v (min/umol)', fontsize=12)
ax2.set_title('Lineweaver-Burk: Inhibition Patterns', fontsize=14)
ax2.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1', fontsize=9)
ax2.tick_params(colors='#cbd5e1')
ax2.xaxis.label.set_color('#cbd5e1')
ax2.yaxis.label.set_color('#cbd5e1')
ax2.title.set_color('#e2e8f0')
for spine in ax2.spines.values():
    spine.set_color('#334155')
ax2.grid(True, alpha=0.15, color='#cbd5e1')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a', edgecolor='none')

print("Enzyme Inhibition Comparison")
print("=" * 45)
print(f"[I] = {I_conc:.1f} mM, Ki = {Ki:.1f} mM")
print(f"alpha = 1 + [I]/Ki = {alpha_c:.2f}")
print()
print(f"{'Type':<20} {'App Km (mM)':<15} {'App Vmax':<15}")
print(f"{'-'*50}")
print(f"{'None':<20} {Km:<15.2f} {Vmax:<15.2f}")
print(f"{'Competitive':<20} {Km*alpha_c:<15.2f} {Vmax:<15.2f}")
print(f"{'Uncompetitive':<20} {Km/alpha_u:<15.2f} {Vmax/alpha_u:<15.2f}")
print(f"{'Mixed (NC)':<20} {Km:<15.2f} {Vmax/alpha_m:<15.2f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Understanding the Results

On the Lineweaver-Burk plot, each inhibition type has a distinct signature: competitive inhibitors change the slope (x-intercept shifts) but share the same y-intercept; uncompetitive inhibitors produce parallel lines (same slope, different intercepts); and mixed/noncompetitive inhibitors change both slope and y-intercept.

Try increasing $[I]$ to see how the curves diverge. At very high substrate concentrations, competitive inhibition can be overcome (curves converge), but uncompetitive and mixed inhibition cannot.

3. Glycolysis Flux Model

Glycolysis converts one molecule of glucose into two molecules of pyruvate through a 10-step enzymatic pathway. The net reaction is:

$$\text{Glucose} + 2\text{NAD}^+ + 2\text{ADP} + 2\text{P}_i \rightarrow 2\text{Pyruvate} + 2\text{NADH} + 2\text{ATP} + 2\text{H}_2\text{O}$$

This simulation models the flux through the three key irreversible regulatory steps: hexokinase (HK), phosphofructokinase-1 (PFK-1), and pyruvate kinase (PK). Each enzyme follows Michaelis-Menten kinetics, and the intermediates glucose-6-phosphate (G6P), fructose-6-phosphate (F6P), fructose-1,6-bisphosphate (FBP), and pyruvate are tracked as a simplified ODE system. PFK-1 is the principal committed step and rate-limiting enzyme, allosterically activated by AMP and inhibited by ATP and citrate.

Python

Parameters

Initial [Glucose] (mM)

Starting glucose concentration

PFK Activity (x)

PFK-1 activity factor (1 = normal)

HK Vmax (x)

Hexokinase Vmax factor

script.py79 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.integrate import odeint

GLUCOSE_VALUE = 5.0
PFK_VALUE = 1.0
HK_VMAX_VALUE = 1.0

glc0 = GLUCOSE_VALUE
pfk_factor = PFK_VALUE
hk_vmax = HK_VMAX_VALUE

# Simplified glycolysis ODE: Glucose -> G6P -> F6P -> FBP -> Pyruvate
# Rate constants (Michaelis-Menten Vmax values, mM/s)
V_hk = 0.8 * hk_vmax     # Hexokinase
K_hk = 0.1               # Km for HK
V_pgi = 2.0              # Phosphoglucose isomerase (fast, near-equilibrium)
K_pgi = 0.4
V_pfk = 0.6 * pfk_factor # PFK-1 (rate-limiting)
K_pfk = 0.2
V_pk = 1.2               # Pyruvate kinase
K_pk = 0.3

def glycolysis(y, t):
    glc, g6p, f6p, fbp, pyr = y
    # Fluxes (Michaelis-Menten)
    J_hk  = V_hk  * glc / (K_hk  + glc)  if glc > 0 else 0
    J_pgi = V_pgi * g6p / (K_pgi + g6p)   if g6p > 0 else 0
    J_pfk = V_pfk * f6p / (K_pfk + f6p)   if f6p > 0 else 0
    J_pk  = V_pk  * fbp / (K_pk  + fbp)   if fbp > 0 else 0

dglc = -J_hk
    dg6p = J_hk - J_pgi
    df6p = J_pgi - J_pfk
    dfbp = J_pfk - J_pk
    dpyr = 2.0 * J_pk  # 2 pyruvate per FBP

return [dglc, dg6p, df6p, dfbp, dpyr]

y0 = [glc0, 0.0, 0.0, 0.0, 0.0]
t = np.linspace(0, 30, 1000)
sol = odeint(glycolysis, y0, t)

fig, ax = plt.subplots(figsize=(10, 7))
fig.patch.set_facecolor('#0f172a')
ax.set_facecolor('#1e293b')

labels = ['Glucose', 'G6P', 'F6P', 'FBP', 'Pyruvate']
colors = ['#4ade80', '#60a5fa', '#fbbf24', '#c084fc', '#f87171']
for i, (lbl, clr) in enumerate(zip(labels, colors)):
    ax.plot(t, sol[:, i], linewidth=2.2, color=clr, label=lbl)

ax.set_xlabel('Time (s)', fontsize=12)
ax.set_ylabel('Concentration (mM)', fontsize=12)
ax.set_title('Glycolysis: Metabolite Concentrations over Time', fontsize=14)
ax.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1', fontsize=10)
ax.tick_params(colors='#cbd5e1')
ax.xaxis.label.set_color('#cbd5e1')
ax.yaxis.label.set_color('#cbd5e1')
ax.title.set_color('#e2e8f0')
for spine in ax.spines.values():
    spine.set_color('#334155')
ax.grid(True, alpha=0.15, color='#cbd5e1')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a', edgecolor='none')

print("Glycolysis Flux Model Results")
print("=" * 40)
print(f"Initial glucose: {glc0:.1f} mM")
print(f"PFK activity factor: {pfk_factor:.2f}")
print(f"HK Vmax factor: {hk_vmax:.2f}")
print(f"Final glucose: {sol[-1,0]:.3f} mM")
print(f"Final pyruvate: {sol[-1,4]:.3f} mM")
print(f"Theoretical max pyruvate: {2*glc0:.1f} mM")
print(f"Conversion efficiency: {sol[-1,4]/(2*glc0)*100:.1f}%")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Understanding the Results

Glucose is consumed by hexokinase and converted through intermediates to pyruvate. Because PFK-1 is the slowest step, you can see F6P accumulates while the PFK activity factor is low. Increasing PFK activity relieves this bottleneck and accelerates the entire pathway.

Note that each glucose yields two pyruvates (the C6 sugar is cleaved into two C3 fragments at the aldolase step). Try lowering PFK activity to 0.2 to simulate citrate inhibition, or raising HK Vmax to simulate insulin-stimulated glucose uptake.

4. TCA Cycle Dynamics

The tricarboxylic acid (TCA) cycle, also known as the citric acid cycle or Krebs cycle, is the central metabolic hub where acetyl-CoA is oxidized to $\text{CO}_2$, generating the reduced cofactors NADH and $\text{FADH}_2$ that drive oxidative phosphorylation. Each turn of the cycle produces:

$$\text{Acetyl-CoA} + 3\text{NAD}^+ + \text{FAD} + \text{GDP} + \text{P}_i \rightarrow 2\text{CO}_2 + 3\text{NADH} + \text{FADH}_2 + \text{GTP} + \text{CoA}$$

This simulation models the concentrations of key TCA intermediates over time as acetyl-CoA enters the cycle. The rate of each step depends on substrate availability and the $\text{NAD}^+/\text{NADH}$ ratio, which controls the oxidative steps catalyzed by isocitrate dehydrogenase, alpha-ketoglutarate dehydrogenase, and malate dehydrogenase.

Python

Parameters

Acetyl-CoA influx (mM/s)

Rate of acetyl-CoA entering the cycle

NAD+/NADH ratio

Oxidized/reduced cofactor ratio

script.py84 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.integrate import odeint

ACOA_INFLUX_VALUE = 0.5
NAD_RATIO_VALUE = 5.0

acoa_influx = ACOA_INFLUX_VALUE
nad_ratio = NAD_RATIO_VALUE  # NAD+/NADH ratio (higher = more oxidized)

# Simplified TCA cycle ODE model
# Species: OAA, citrate, isocitrate, alpha-KG, succinate, malate
# Rate constants scaled by NAD+/NADH ratio for oxidative steps

k_cs  = 1.0    # citrate synthase
K_cs  = 0.2
k_aco = 3.0    # aconitase (fast equilibrium)
K_aco = 0.3
k_idh = 0.8 * (nad_ratio / 5.0)  # isocitrate dehydrogenase (NAD-dependent)
K_idh = 0.15
k_kgdh = 0.6 * (nad_ratio / 5.0) # alpha-KG dehydrogenase (NAD-dependent)
K_kgdh = 0.2
k_sdh = 1.5    # succinate dehydrogenase (FAD)
K_sdh = 0.3
k_mdh = 0.7 * (nad_ratio / 5.0)  # malate dehydrogenase (NAD-dependent)
K_mdh = 0.25

def tca_cycle(y, t):
    oaa, cit, iso, akg, suc, mal = y

J_cs   = k_cs  * oaa / (K_cs + oaa)  * (acoa_influx / (0.1 + acoa_influx))
    J_aco  = k_aco * cit / (K_aco + cit) if cit > 0 else 0
    J_idh  = k_idh * iso / (K_idh + iso) if iso > 0 else 0
    J_kgdh = k_kgdh * akg / (K_kgdh + akg) if akg > 0 else 0
    J_sdh  = k_sdh * suc / (K_sdh + suc) if suc > 0 else 0
    J_mdh  = k_mdh * mal / (K_mdh + mal) if mal > 0 else 0

doaa = J_mdh - J_cs
    dcit = J_cs - J_aco
    diso = J_aco - J_idh
    dakg = J_idh - J_kgdh
    dsuc = J_kgdh - J_sdh
    dmal = J_sdh - J_mdh

return [doaa, dcit, diso, dakg, dsuc, dmal]

y0 = [0.1, 0.05, 0.02, 0.02, 0.05, 0.05]
t = np.linspace(0, 40, 1000)
sol = odeint(tca_cycle, y0, t)

fig, ax = plt.subplots(figsize=(10, 7))
fig.patch.set_facecolor('#0f172a')
ax.set_facecolor('#1e293b')

names = ['OAA', 'Citrate', 'Isocitrate', 'alpha-KG', 'Succinate', 'Malate']
colors = ['#f87171', '#4ade80', '#60a5fa', '#fbbf24', '#c084fc', '#fb923c']
for i, (nm, clr) in enumerate(zip(names, colors)):
    ax.plot(t, sol[:, i], linewidth=2.2, color=clr, label=nm)

ax.set_xlabel('Time (s)', fontsize=12)
ax.set_ylabel('Concentration (mM)', fontsize=12)
ax.set_title('TCA Cycle Intermediate Dynamics', fontsize=14)
ax.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1', fontsize=10, ncol=2)
ax.tick_params(colors='#cbd5e1')
ax.xaxis.label.set_color('#cbd5e1')
ax.yaxis.label.set_color('#cbd5e1')
ax.title.set_color('#e2e8f0')
for spine in ax.spines.values():
    spine.set_color('#334155')
ax.grid(True, alpha=0.15, color='#cbd5e1')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a', edgecolor='none')

print("TCA Cycle Dynamics")
print("=" * 40)
print(f"Acetyl-CoA influx: {acoa_influx:.2f} mM/s")
print(f"NAD+/NADH ratio: {nad_ratio:.1f}")
print(f"\nSteady-state concentrations (approx at t={t[-1]:.0f}s):")
for nm, val in zip(names, sol[-1]):
    print(f"  {nm:<15} {val:.4f} mM")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Understanding the Results

The cycle reaches a quasi-steady state where intermediate concentrations stabilize. With a high $\text{NAD}^+/\text{NADH}$ ratio, the three NAD-dependent dehydrogenases are accelerated, draining isocitrate, alpha-KG, and malate faster. Lowering this ratio (simulating reduced oxygen supply) causes these intermediates to accumulate.

OAA concentration is critical because it must combine with acetyl-CoA in the citrate synthase step. If OAA is depleted (anaplerotic reactions are insufficient), the cycle stalls regardless of acetyl-CoA availability.

5. Electron Transport Chain & ATP Yield Calculator

The electron transport chain (ETC) couples the oxidation of NADH and $\text{FADH}_2$ to the generation of a proton gradient across the inner mitochondrial membrane. The free energy released by electron transfer through Complexes I, III, and IV is:

$$\Delta G = -nF\Delta E'_0, \quad \text{where } \Delta E'_0(\text{NADH} \to \text{O}_2) = 1.14\text{ V}$$

The P/O ratio defines how many ATP are produced per oxygen atom reduced. For NADH (entering at Complex I), the theoretical P/O ratio is approximately 2.5; for $\text{FADH}_2$ (entering at Complex II), it is approximately 1.5. This Fortran program calculates total ATP yield and thermodynamic efficiency based on the standard free energy of glucose oxidation ($\Delta G^\circ = -2870$ kJ/mol) and ATP hydrolysis ($\Delta G^\circ = -30.5$ kJ/mol).

Fortran

Parameters

NADH

Number of NADH molecules

FADH2

Number of FADH2 molecules

P/O (NADH)

ATP per 1/2 O2 from NADH

P/O (FADH2)

ATP per 1/2 O2 from FADH2

Substrate ATP

Substrate-level phosphorylation ATP

program.f9053 lines

program atp_yield
  implicit none
  real :: nadh_val, fadh2_val, po_nadh, po_fadh2, substrate_atp
  real :: atp_from_nadh, atp_from_fadh2, total_atp, efficiency
  real :: delta_g_glucose, delta_g_atp

nadh_val = NADH_VALUE
  fadh2_val = FADH2_VALUE
  po_nadh = PO_NADH_VALUE
  po_fadh2 = PO_FADH2_VALUE
  substrate_atp = SUBSTRATE_ATP_VALUE

delta_g_glucose = 2870.0
  delta_g_atp = 30.5

atp_from_nadh = nadh_val * po_nadh
  atp_from_fadh2 = fadh2_val * po_fadh2
  total_atp = atp_from_nadh + atp_from_fadh2 + substrate_atp
  efficiency = (total_atp * delta_g_atp) / delta_g_glucose * 100.0

write(*,'(A)') '================================================'
  write(*,'(A)') '    Electron Transport Chain & ATP Yield'
  write(*,'(A)') '================================================'
  write(*,'(A)') ''
  write(*,'(A)') '--- Inputs ---'
  write(*,'(A,F8.1)') 'NADH molecules:          ', nadh_val
  write(*,'(A,F8.1)') 'FADH2 molecules:         ', fadh2_val
  write(*,'(A,F8.2)') 'P/O ratio (NADH):        ', po_nadh
  write(*,'(A,F8.2)') 'P/O ratio (FADH2):       ', po_fadh2
  write(*,'(A,F8.1)') 'Substrate-level ATP:     ', substrate_atp
  write(*,'(A)') ''
  write(*,'(A)') '--- ATP Breakdown ---'
  write(*,'(A,F8.1)') 'ATP from NADH:           ', atp_from_nadh
  write(*,'(A,F8.1)') 'ATP from FADH2:          ', atp_from_fadh2
  write(*,'(A,F8.1)') 'Substrate-level ATP:     ', substrate_atp
  write(*,'(A)') '                         --------'
  write(*,'(A,F8.1)') 'TOTAL ATP yield:         ', total_atp
  write(*,'(A)') ''
  write(*,'(A)') '--- Thermodynamics ---'
  write(*,'(A,F8.1,A)') 'Energy captured:         ', &
       total_atp * delta_g_atp, ' kJ/mol'
  write(*,'(A,F8.1,A)') 'Glucose free energy:     ', &
       delta_g_glucose, ' kJ/mol'
  write(*,'(A,F8.1,A)') 'Thermodynamic efficiency:', efficiency, ' %'
  write(*,'(A)') ''
  write(*,'(A)') '--- Per Substrate Comparison ---'
  write(*,'(A)') 'Glucose (complete oxidation): 10 NADH + 2 FADH2 + 4 ATP'
  write(*,'(A,F6.1,A)') '  -> ', 10.0*po_nadh + 2.0*po_fadh2 + 4.0, ' ATP'
  write(*,'(A)') 'Palmitate (16C fatty acid):   ~31 NADH + 15 FADH2 + 1 ATP (net)'
  write(*,'(A,F6.1,A)') '  -> ', 31.0*po_nadh + 15.0*po_fadh2 + 1.0, ' ATP'
  write(*,'(A)') '================================================'

end program atp_yield

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

Understanding the Results

For complete glucose oxidation (glycolysis + pyruvate dehydrogenase + TCA cycle), the default values yield 10 NADH, 2 $\text{FADH}_2$, and 4 substrate-level ATP (including 2 GTP from succinyl-CoA synthetase). With standard P/O ratios of 2.5 and 1.5, this gives approximately 30-32 ATP per glucose.

The thermodynamic efficiency of about 34% reflects the fraction of glucose free energy captured as ATP. The remainder is released as heat, which is essential for maintaining body temperature. Try comparing glucose vs palmitate by adjusting to 31 NADH, 15 $\text{FADH}_2$, and 1 substrate-level ATP.

6. Beta-Oxidation Energy Calculator

Beta-oxidation is the catabolic process by which fatty acyl-CoA is sequentially shortened by two carbons per cycle, producing acetyl-CoA, NADH, and $\text{FADH}_2$ at each step. For a saturated fatty acid with $n$ carbon atoms, the number of beta-oxidation cycles is $n/2 - 1$, yielding:

$$\text{Acetyl-CoA} = \frac{n}{2}, \quad \text{NADH} = \frac{n}{2}-1, \quad \text{FADH}_2 = \frac{n}{2}-1$$

Each acetyl-CoA entering the TCA cycle yields 3 NADH + 1 $\text{FADH}_2$ + 1 GTP, so the total ATP yield from a single fatty acid molecule is substantial. However, activation to fatty acyl-CoA costs 2 ATP equivalents. This program calculates the complete energy balance and displays a breakdown by source.

Python

Parameters

Chain length (carbons)

Even-numbered carbon chain length

script.py114 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

CHAIN_LENGTH_VALUE = 16

n = int(CHAIN_LENGTH_VALUE)
if n % 2 != 0:
    n = n - 1  # ensure even chain length

cycles = n // 2 - 1
acetyl_coa = n // 2

# Beta-oxidation products per cycle
nadh_beta = cycles
fadh2_beta = cycles

# TCA cycle products per acetyl-CoA
nadh_tca = 3 * acetyl_coa
fadh2_tca = 1 * acetyl_coa
gtp_tca = 1 * acetyl_coa

# Total reduced cofactors
total_nadh = nadh_beta + nadh_tca
total_fadh2 = fadh2_beta + fadh2_tca

# ATP yield (P/O = 2.5 for NADH, 1.5 for FADH2)
atp_nadh = total_nadh * 2.5
atp_fadh2 = total_fadh2 * 1.5
atp_gtp = gtp_tca * 1.0
activation_cost = 2.0  # 2 ATP equivalents for activation

total_atp = atp_nadh + atp_fadh2 + atp_gtp - activation_cost

# Bar chart breakdown
categories = ['NADH\n(beta-ox)', 'FADH2\n(beta-ox)', 'NADH\n(TCA)', 'FADH2\n(TCA)', 'GTP\n(TCA)', 'Activation\ncost']
atp_values = [nadh_beta * 2.5, fadh2_beta * 1.5, nadh_tca * 2.5, fadh2_tca * 1.5, gtp_tca, -activation_cost]
colors = ['#4ade80', '#60a5fa', '#22d3ee', '#a78bfa', '#fbbf24', '#f87171']

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
fig.patch.set_facecolor('#0f172a')

# Left: Bar chart
ax1.set_facecolor('#1e293b')
bars = ax1.bar(categories, atp_values, color=colors, edgecolor='#334155', linewidth=1.2)
for bar, val in zip(bars, atp_values):
    ax1.text(bar.get_x() + bar.get_width()/2, bar.get_height() + (1 if val >= 0 else -3),
             f'{val:+.1f}', ha='center', va='bottom' if val >= 0 else 'top',
             color='#e2e8f0', fontsize=10, fontweight='bold')
ax1.axhline(y=0, color='#475569', linewidth=0.8)
ax1.set_ylabel('ATP equivalents', fontsize=12)
ax1.set_title(f'ATP Breakdown: C{n} Fatty Acid', fontsize=14)
ax1.tick_params(colors='#cbd5e1', labelsize=9)
ax1.xaxis.label.set_color('#cbd5e1')
ax1.yaxis.label.set_color('#cbd5e1')
ax1.title.set_color('#e2e8f0')
for spine in ax1.spines.values():
    spine.set_color('#334155')

# Right: Comparison across chain lengths
ax2.set_facecolor('#1e293b')
chain_lengths = np.arange(4, 24, 2)
net_atp_list = []
atp_per_carbon = []
for c in chain_lengths:
    cy = c // 2 - 1
    ac = c // 2
    total_n = (cy + 3*ac) * 2.5 + (cy + ac) * 1.5 + ac - 2
    net_atp_list.append(total_n)
    atp_per_carbon.append(total_n / c)

ax2.bar(chain_lengths - 0.4, net_atp_list, width=0.8, color='#4ade80', alpha=0.8, label='Total ATP')
ax2_twin = ax2.twinx()
ax2_twin.plot(chain_lengths, atp_per_carbon, 'o-', color='#f87171', linewidth=2, markersize=6, label='ATP/carbon')
ax2_twin.set_ylabel('ATP per carbon', fontsize=12)
ax2_twin.yaxis.label.set_color('#f87171')
ax2_twin.tick_params(colors='#f87171')

ax2.set_xlabel('Chain length (carbons)', fontsize=12)
ax2.set_ylabel('Net ATP yield', fontsize=12)
ax2.set_title('ATP Yield vs Chain Length', fontsize=14)
ax2.tick_params(colors='#cbd5e1')
ax2.xaxis.label.set_color('#cbd5e1')
ax2.yaxis.label.set_color('#cbd5e1')
ax2.title.set_color('#e2e8f0')
for spine in ax2.spines.values():
    spine.set_color('#334155')
for spine in ax2_twin.spines.values():
    spine.set_color('#334155')

lines1, labels1 = ax2.get_legend_handles_labels()
lines2, labels2 = ax2_twin.get_legend_handles_labels()
ax2.legend(lines1+lines2, labels1+labels2, facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1', fontsize=9)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a', edgecolor='none')

print(f"Beta-Oxidation Energy Balance: C{n} Fatty Acid")
print("=" * 50)
print(f"Carbon chain length: {n}")
print(f"Beta-oxidation cycles: {cycles}")
print(f"Acetyl-CoA produced: {acetyl_coa}")
print(f"")
print(f"NADH (beta-ox): {nadh_beta}  -> {nadh_beta*2.5:.1f} ATP")
print(f"FADH2 (beta-ox): {fadh2_beta}  -> {fadh2_beta*1.5:.1f} ATP")
print(f"NADH (TCA): {nadh_tca}  -> {nadh_tca*2.5:.1f} ATP")
print(f"FADH2 (TCA): {fadh2_tca}  -> {fadh2_tca*1.5:.1f} ATP")
print(f"GTP (TCA): {gtp_tca}  -> {gtp_tca:.1f} ATP")
print(f"Activation cost: -2.0 ATP")
print(f"")
print(f"NET ATP YIELD: {total_atp:.1f}")
print(f"ATP per carbon: {total_atp/n:.2f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Understanding the Results

Palmitate (C16) yields approximately 106 net ATP, far exceeding the ~30 ATP from glucose. This is why fats are the most energy-dense macronutrient (9 kcal/g vs 4 kcal/g for carbohydrates). The right panel shows that ATP per carbon increases with chain length because the fixed 2 ATP activation cost becomes proportionally smaller.

The majority of ATP comes from the TCA cycle NADH, reflecting the fact that complete oxidation of acetyl-CoA through the cycle and ETC is the dominant energy source.

7. Urea Cycle & Nitrogen Balance

The urea cycle is the primary pathway for disposing of excess nitrogen from amino acid catabolism. It converts toxic ammonia ($\text{NH}_3$) and the amino group of aspartate into urea, which is excreted by the kidneys. The net reaction is:

$$\text{NH}_3 + \text{CO}_2 + \text{Aspartate} + 3\text{ATP} \rightarrow \text{Urea} + \text{Fumarate} + 2\text{ADP} + \text{AMP} + 2\text{P}_i + \text{PP}_i$$

This Fortran program calculates nitrogen balance based on dietary protein intake. Since proteins average about 16% nitrogen by mass, and each gram of nitrogen produces approximately$2.14$ g of urea (MW urea = 60, MW N = 14, ratio = 60/28), we can estimate daily urea production and the ATP cost of nitrogen disposal.

Fortran

Parameters

Protein intake (g/day)

Daily dietary protein

Oxidation fraction

Fraction of amino acids oxidized (vs used for synthesis)

program.f9071 lines

program urea_cycle
  implicit none
  real :: protein_intake, oxidation_frac
  real :: nitrogen_in, nitrogen_oxidized, urea_produced
  real :: ammonia_load, atp_cost, urea_nitrogen
  real :: protein_for_synthesis, nitrogen_retained

protein_intake = PROTEIN_VALUE
  oxidation_frac = OXIDATION_VALUE

! Nitrogen calculations
  nitrogen_in = protein_intake * 0.16
  nitrogen_oxidized = nitrogen_in * oxidation_frac
  nitrogen_retained = nitrogen_in * (1.0 - oxidation_frac)
  protein_for_synthesis = nitrogen_retained / 0.16

! Each urea molecule contains 2 nitrogen atoms
  ! MW urea = 60 g/mol, MW N = 14 g/mol
  ! 1 g N -> (60 / 28) = 2.143 g urea
  urea_produced = nitrogen_oxidized * (60.0 / 28.0)
  urea_nitrogen = nitrogen_oxidized

! Ammonia load (intermediate)
  ammonia_load = nitrogen_oxidized * (17.0 / 14.0)

! ATP cost: 4 high-energy phosphate bonds per urea
  ! (3 ATP -> 2 ADP + AMP, but AMP->ADP costs 1 ATP equivalent)
  ! Moles of N oxidized
  ! nitrogen_oxidized is in grams, MW N = 14
  atp_cost = (nitrogen_oxidized / 14.0) * 0.5 * 4.0
  ! (0.5 because 2 N per urea, times 4 ATP per urea)

write(*,'(A)') '================================================'
  write(*,'(A)') '     Urea Cycle & Nitrogen Balance Calculator'
  write(*,'(A)') '================================================'
  write(*,'(A)') ''
  write(*,'(A)') '--- Dietary Input ---'
  write(*,'(A,F8.1,A)') 'Protein intake:       ', protein_intake, ' g/day'
  write(*,'(A,F8.1,A)') 'Oxidation fraction:   ', oxidation_frac*100, ' %'
  write(*,'(A)') ''
  write(*,'(A)') '--- Nitrogen Balance ---'
  write(*,'(A,F8.2,A)') 'Total nitrogen in:    ', nitrogen_in, ' g/day'
  write(*,'(A,F8.2,A)') 'Nitrogen oxidized:    ', nitrogen_oxidized, ' g/day'
  write(*,'(A,F8.2,A)') 'Nitrogen retained:    ', nitrogen_retained, ' g/day'
  write(*,'(A,F8.1,A)') 'Protein synthesis:    ', protein_for_synthesis, ' g/day'
  write(*,'(A)') ''
  write(*,'(A)') '--- Urea Production ---'
  write(*,'(A,F8.2,A)') 'Urea produced:        ', urea_produced, ' g/day'
  write(*,'(A,F8.2,A)') 'Urea nitrogen:        ', urea_nitrogen, ' g/day'
  write(*,'(A,F8.2,A)') 'Ammonia intermediate: ', ammonia_load, ' g/day'
  write(*,'(A)') ''
  write(*,'(A)') '--- Energy Cost ---'
  write(*,'(A,F8.2,A)') 'Urea (moles/day):     ', &
       urea_produced / 60.0, ' mol/day'
  write(*,'(A,F8.1,A)') 'ATP cost:             ', atp_cost, ' mol/day'
  write(*,'(A,F8.1,A)') 'Energy cost (~30.5 kJ/mol ATP): ', &
       atp_cost * 30.5, ' kJ/day'
  write(*,'(A)') ''
  write(*,'(A)') '--- Clinical Notes ---'
  if (nitrogen_oxidized .gt. 16.0) then
    write(*,'(A)') 'WARNING: Very high nitrogen load.'
    write(*,'(A)') 'Risk of hyperammonemia if urea cycle impaired.'
  else if (nitrogen_oxidized .lt. 3.0) then
    write(*,'(A)') 'Low nitrogen oxidation - possible protein-sparing'
    write(*,'(A)') 'or inadequate protein intake.'
  else
    write(*,'(A)') 'Nitrogen load within normal range.'
  end if
  write(*,'(A)') '================================================'

end program urea_cycle

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

Understanding the Results

On a typical 80 g/day protein diet with 70% oxidation, the body must dispose of about 9 g of nitrogen daily as urea. This costs approximately 4 ATP per urea molecule, representing a non-trivial metabolic expense of nitrogen disposal.

In urea cycle disorders (e.g., ornithine transcarbamylase deficiency), ammonia accumulates to toxic levels, causing encephalopathy. Treatment involves reducing protein intake and providing alternative nitrogen disposal pathways (sodium benzoate, phenylbutyrate). Try increasing protein intake to see how the nitrogen load scales.

8. Blood Glucose Homeostasis

Blood glucose is maintained within a narrow range (4-6 mM fasting) by the opposing actions of insulin and glucagon. After a meal, rising glucose stimulates pancreatic beta-cells to secrete insulin, which promotes glucose uptake into muscle and adipose tissue and stimulates glycogen synthesis. Between meals, falling glucose triggers alpha-cell glucagon secretion, which activates hepatic glycogenolysis and gluconeogenesis.

This model uses a system of ODEs to simulate the glucose-insulin-glucagon dynamics following a glucose load. In type 2 diabetes, insulin resistance reduces the effectiveness of insulin, requiring higher insulin levels to achieve the same glucose disposal. The model captures this by scaling the insulin sensitivity parameter $\sigma$:

$$\frac{d[\text{Glucose}]}{dt} = R_{\text{meal}}(t) + R_{\text{glucagon}} - \sigma \cdot [\text{Insulin}] \cdot [\text{Glucose}] - k_{\text{basal}} \cdot [\text{Glucose}]$$

Python

Parameters

Glucose load (mg)

Oral glucose load (75g = standard OGTT)

Insulin sensitivity

1.0 = normal, <1 = resistant

Basal glucagon

Relative basal glucagon level

script.py94 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.integrate import odeint

MEAL_VALUE = 75.0
SENSITIVITY_VALUE = 1.0
GLUCAGON_VALUE = 1.0

meal_glucose = MEAL_VALUE   # mg glucose load
sensitivity = SENSITIVITY_VALUE  # insulin sensitivity (1 = normal, <1 = resistant)
glucagon_basal = GLUCAGON_VALUE   # relative basal glucagon

def glucose_model(y, t, sigma):
    G, I, Gn = y  # glucose (mM), insulin (uU/mL), glucagon (pg/mL)

# Meal absorption: first-order absorption from gut
    R_meal = (meal_glucose / 180.0 * 1000 / 5000) * 2.0 * np.exp(-0.03 * t) * (0.03 * t) if t < 300 else 0

# Insulin secretion: sigmoid response to glucose
    I_secretion = 50.0 * (G**2) / (5.0**2 + G**2) - 0.2 * I

# Glucagon: suppressed by glucose and insulin
    Gn_secretion = glucagon_basal * 40.0 / (1.0 + G / 4.0) - 0.15 * Gn

# Glucose dynamics
    R_glucagon = 0.005 * Gn  # hepatic glucose output
    disposal = sigma * 0.004 * I * G  # insulin-mediated uptake
    basal_uptake = 0.01 * G  # insulin-independent uptake

dG = R_meal + R_glucagon - disposal - basal_uptake
    dI = I_secretion
    dGn = Gn_secretion

return [dG, dI, dGn]

t = np.linspace(0, 360, 2000)  # 6 hours in minutes

# Normal case
y0_normal = [5.0, 10.0, 50.0]
sol_normal = odeint(glucose_model, y0_normal, t, args=(sensitivity,))

# Type 2 diabetes (reduced sensitivity)
sol_t2d = odeint(glucose_model, y0_normal, t, args=(sensitivity * 0.3,))

fig, axes = plt.subplots(3, 1, figsize=(10, 10))
fig.patch.set_facecolor('#0f172a')

titles = ['Blood Glucose (mM)', 'Plasma Insulin (uU/mL)', 'Plasma Glucagon (pg/mL)']
colors_n = ['#4ade80', '#60a5fa', '#fbbf24']
colors_d = ['#f87171', '#c084fc', '#fb923c']
ylabels = ['Glucose (mM)', 'Insulin (uU/mL)', 'Glucagon (pg/mL)']

for i, ax in enumerate(axes):
    ax.set_facecolor('#1e293b')
    ax.plot(t, sol_normal[:, i], linewidth=2.2, color=colors_n[i], label='Normal')
    ax.plot(t, sol_t2d[:, i], linewidth=2.2, color=colors_d[i], linestyle='--', label='Type 2 DM')
    ax.set_ylabel(ylabels[i], fontsize=11)
    ax.set_title(titles[i], fontsize=13)
    ax.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1', fontsize=9)
    ax.tick_params(colors='#cbd5e1')
    ax.xaxis.label.set_color('#cbd5e1')
    ax.yaxis.label.set_color('#cbd5e1')
    ax.title.set_color('#e2e8f0')
    for spine in ax.spines.values():
        spine.set_color('#334155')
    ax.grid(True, alpha=0.15, color='#cbd5e1')

# Add reference ranges
axes[0].axhline(y=7.0, color='#f87171', alpha=0.3, linestyle=':', linewidth=1)
axes[0].axhline(y=4.0, color='#fbbf24', alpha=0.3, linestyle=':', linewidth=1)
axes[0].text(350, 7.2, 'Hyperglycemia', color='#f87171', fontsize=8, ha='right', alpha=0.7)
axes[0].text(350, 3.5, 'Hypoglycemia', color='#fbbf24', fontsize=8, ha='right', alpha=0.7)

axes[2].set_xlabel('Time (min)', fontsize=12)
axes[2].xaxis.label.set_color('#cbd5e1')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a', edgecolor='none')

print("Blood Glucose Homeostasis Model")
print("=" * 45)
print(f"Glucose load: {meal_glucose:.0f} mg")
print(f"Insulin sensitivity: {sensitivity:.2f}")
print(f"Basal glucagon factor: {glucagon_basal:.2f}")
print()
print(f"{'Parameter':<25} {'Normal':<15} {'Type 2 DM':<15}")
print(f"{'-'*55}")
print(f"{'Peak glucose (mM)':<25} {max(sol_normal[:,0]):<15.2f} {max(sol_t2d[:,0]):<15.2f}")
print(f"{'Peak insulin (uU/mL)':<25} {max(sol_normal[:,1]):<15.2f} {max(sol_t2d[:,1]):<15.2f}")
print(f"{'Time to peak glc (min)':<25} {t[np.argmax(sol_normal[:,0])]:<15.1f} {t[np.argmax(sol_t2d[:,0])]:<15.1f}")
print(f"{'Glucose at 120 min (mM)':<25} {sol_normal[int(120/360*2000),0]:<15.2f} {sol_t2d[int(120/360*2000),0]:<15.2f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Understanding the Results

In the normal response (solid lines), glucose rises after the meal, peaks around 30-60 minutes, then returns to baseline within 2 hours as insulin-mediated uptake clears it. Insulin mirrors the glucose curve with a slight delay, while glucagon is suppressed during the postprandial period.

In the type 2 diabetes simulation (dashed lines, 30% of normal sensitivity), glucose rises higher and takes much longer to return to baseline. The beta-cells compensate by secreting more insulin (hyperinsulinemia), but the reduced tissue response means glucose disposal is impaired. This is the hallmark of insulin resistance seen in the early stages of type 2 diabetes.

9. Metabolic Flux Analysis

Metabolic flux analysis (MFA) determines the rates of metabolic reactions in a network at steady state. At any metabolic branch point, the flux must be conserved (Kirchhoff-like law). For glucose entering a cell, the three major fates are glycolysis (energy production), the pentose phosphate pathway (PPP, for NADPH and ribose-5-phosphate), and glycogenesis (storage). The fluxes must satisfy:

$$J_{\text{uptake}} = J_{\text{glycolysis}} + J_{\text{PPP}} + J_{\text{glycogenesis}}$$

This simulation solves a simple linear program to find the optimal flux distribution given constraints on ATP demand and NADPH demand. Glycolysis yields 2 ATP per glucose; the PPP yields 2 NADPH per glucose-6-phosphate entering the oxidative branch; glycogenesis consumes 1 ATP per glucose stored. The solver maximizes glucose storage while meeting energy and biosynthetic demands.

Python

Parameters

Glucose uptake (mmol/h)

Total glucose uptake rate

ATP demand (mmol/h)

Cellular ATP requirement

NADPH demand (mmol/h)

NADPH requirement (biosynthesis)

script.py124 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

UPTAKE_VALUE = 10.0
ATP_DEMAND_VALUE = 12.0
NADPH_DEMAND_VALUE = 4.0

J_uptake = UPTAKE_VALUE    # glucose uptake (mmol/h)
ATP_demand = ATP_DEMAND_VALUE   # ATP demand (mmol/h)
NADPH_demand = NADPH_DEMAND_VALUE  # NADPH demand (mmol/h)

# Yields per glucose unit
# Glycolysis: 2 ATP, 0 NADPH
# PPP: 0 net ATP, 2 NADPH (oxidative branch)
# Glycogenesis: -1 ATP (costs 1 UTP = 1 ATP equivalent), 0 NADPH

# Solve: minimize waste / maximize glycogenesis
# J_glyc + J_ppp + J_glycogen = J_uptake
# 2*J_glyc - 1*J_glycogen >= ATP_demand
# 2*J_ppp >= NADPH_demand

# Minimum PPP flux
J_ppp = max(NADPH_demand / 2.0, 0)

# Remaining after PPP
remaining = J_uptake - J_ppp

# Need: 2*J_glyc - J_glycogen >= ATP_demand, and J_glyc + J_glycogen = remaining
# J_glyc = remaining - J_glycogen
# 2*(remaining - J_glycogen) - J_glycogen >= ATP_demand
# 2*remaining - 3*J_glycogen >= ATP_demand
# J_glycogen <= (2*remaining - ATP_demand) / 3

J_glycogen_max = (2 * remaining - ATP_demand) / 3.0
J_glycogen = max(J_glycogen_max, 0)
J_glyc = remaining - J_glycogen

# Check if feasible
atp_produced = 2 * J_glyc - J_glycogen
nadph_produced = 2 * J_ppp
feasible = atp_produced >= ATP_demand * 0.99 and J_glyc >= 0 and J_ppp >= 0 and J_glycogen >= 0

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
fig.patch.set_facecolor('#0f172a')

# Left: Pie chart of flux distribution
ax1.set_facecolor('#0f172a')
sizes = [J_glyc, J_ppp, J_glycogen]
labels_pie = [f'Glycolysis\n{J_glyc:.2f}', f'PPP\n{J_ppp:.2f}', f'Glycogenesis\n{J_glycogen:.2f}']
colors_pie = ['#4ade80', '#60a5fa', '#fbbf24']
explode = (0.05, 0.05, 0.05)

# Filter out zero or negative values
valid = [(s, l, c, e) for s, l, c, e in zip(sizes, labels_pie, colors_pie, explode) if s > 0.01]
if valid:
    v_sizes, v_labels, v_colors, v_explode = zip(*valid)
    wedges, texts, autotexts = ax1.pie(v_sizes, labels=v_labels, colors=v_colors,
                                         explode=v_explode, autopct='%1.1f%%',
                                         textprops={'color': '#e2e8f0', 'fontsize': 11},
                                         pctdistance=0.6)
    for at in autotexts:
        at.set_color('#1e293b')
        at.set_fontweight('bold')
else:
    ax1.text(0.5, 0.5, 'Infeasible', transform=ax1.transAxes, ha='center', color='#f87171', fontsize=18)

ax1.set_title(f'Flux Distribution (J_uptake = {J_uptake:.1f} mmol/h)', fontsize=13, color='#e2e8f0')

# Right: Sensitivity bar chart - vary ATP demand
ax2.set_facecolor('#1e293b')
atp_range = np.linspace(0, 2 * J_uptake, 50)
j_glyc_arr = []
j_ppp_arr = []
j_glycogen_arr = []

for atp_d in atp_range:
    jp = max(NADPH_demand / 2.0, 0)
    rem = J_uptake - jp
    jg_max = (2 * rem - atp_d) / 3.0
    jg = max(jg_max, 0)
    jgl = rem - jg
    if jgl < 0:
        jgl = rem
        jg = 0
    j_glyc_arr.append(jgl)
    j_ppp_arr.append(jp)
    j_glycogen_arr.append(jg)

ax2.stackplot(atp_range, j_glyc_arr, j_ppp_arr, j_glycogen_arr,
              colors=['#4ade80', '#60a5fa', '#fbbf24'], alpha=0.8,
              labels=['Glycolysis', 'PPP', 'Glycogenesis'])
ax2.axvline(x=ATP_demand, color='#f87171', linestyle='--', linewidth=1.5, label=f'Current ATP demand')
ax2.set_xlabel('ATP demand (mmol/h)', fontsize=12)
ax2.set_ylabel('Flux (mmol glucose/h)', fontsize=12)
ax2.set_title('Flux Distribution vs ATP Demand', fontsize=13)
ax2.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1', fontsize=9, loc='upper right')
ax2.tick_params(colors='#cbd5e1')
ax2.xaxis.label.set_color('#cbd5e1')
ax2.yaxis.label.set_color('#cbd5e1')
ax2.title.set_color('#e2e8f0')
for spine in ax2.spines.values():
    spine.set_color('#334155')
ax2.grid(True, alpha=0.15, color='#cbd5e1')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a', edgecolor='none')

print("Metabolic Flux Analysis")
print("=" * 45)
print(f"Glucose uptake: {J_uptake:.2f} mmol/h")
print(f"ATP demand: {ATP_demand:.2f} mmol/h")
print(f"NADPH demand: {NADPH_demand:.2f} mmol/h")
print()
print("Optimal flux distribution:")
print(f"  Glycolysis:    {J_glyc:.3f} mmol/h ({J_glyc/J_uptake*100:.1f}%)")
print(f"  PPP:           {J_ppp:.3f} mmol/h ({J_ppp/J_uptake*100:.1f}%)")
print(f"  Glycogenesis:  {J_glycogen:.3f} mmol/h ({J_glycogen/J_uptake*100:.1f}%)")
print()
print(f"ATP produced: {atp_produced:.2f} mmol/h (demand: {ATP_demand:.2f})")
print(f"NADPH produced: {nadph_produced:.2f} mmol/h (demand: {NADPH_demand:.2f})")
print(f"Feasible: {'Yes' if feasible else 'No - demands exceed capacity'}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Understanding the Results

The pie chart shows how glucose is partitioned among the three pathways. The solver first allocates enough flux to the PPP to meet NADPH demand, then distributes the remainder between glycolysis and glycogenesis to meet ATP demand while maximizing storage.

The right panel shows how the flux distribution shifts as ATP demand increases: glycolysis claims a larger share at the expense of glycogenesis. When ATP demand exceeds the capacity of available glucose, glycogen storage drops to zero and all glucose flows through glycolysis. This mirrors the fed-to-fasting transition in vivo.

10. The Warburg Effect

Otto Warburg observed in 1924 that cancer cells preferentially use glycolysis for ATP production even in the presence of oxygen — a phenomenon called aerobic glycolysis or the Warburg effect. While oxidative phosphorylation yields ~36 ATP per glucose, aerobic glycolysis yields only 2 ATP but proceeds much faster and provides biosynthetic precursors.

The tradeoff can be understood quantitatively. If the glycolytic rate is $r_{\text{glyc}}$ (glucose/s) and oxidative phosphorylation rate is $r_{\text{oxphos}}$, then the total ATP production rate for a cell using fraction $f$ of glucose glycolytically is:

$$R_{\text{ATP}} = 2 \cdot f \cdot J_{\text{glc}} + 36 \cdot (1 - f) \cdot J_{\text{glc}} \cdot \frac{[\text{O}_2]}{K_{\text{O}_2} + [\text{O}_2]}$$

This simulation compares the metabolic strategies, showing that while glycolysis is less efficient per glucose, it can produce ATP faster when glucose is abundant. Cancer cells exploit this by upregulating glucose transporters (GLUT1) and glycolytic enzymes.

Python

Parameters

Glucose uptake (a.u.)

Glucose uptake rate

O2 availability

Relative oxygen availability (1 = normoxia)

Glycolytic fraction

Fraction of glucose going to glycolysis

script.py142 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

GLC_UPTAKE_VALUE = 10.0
O2_AVAIL_VALUE = 1.0
F_GLYC_VALUE = 0.5

J_glc = GLC_UPTAKE_VALUE      # glucose uptake rate (arbitrary units)
O2_availability = O2_AVAIL_VALUE  # relative O2 (1 = normoxia)
f_glyc_default = F_GLYC_VALUE  # default glycolytic fraction

# O2-dependent oxphos efficiency
K_O2 = 0.2  # half-saturation for O2

f_range = np.linspace(0, 1, 500)

# ATP production rate
atp_glyc = 2.0 * f_range * J_glc
atp_oxphos = 36.0 * (1 - f_range) * J_glc * O2_availability / (K_O2 + O2_availability)
atp_total = atp_glyc + atp_oxphos

# Lactate production (2 per glucose going through glycolysis)
lactate = 2.0 * f_range * J_glc

# Biomass precursor production (glycolytic intermediates diverted)
biomass_score = f_range * J_glc * 0.5  # proportional to glycolytic flux

fig, axes = plt.subplots(2, 2, figsize=(12, 10))
fig.patch.set_facecolor('#0f172a')

# Top-left: ATP production breakdown
ax = axes[0, 0]
ax.set_facecolor('#1e293b')
ax.fill_between(f_range, 0, atp_oxphos, color='#60a5fa', alpha=0.6, label='OxPhos ATP')
ax.fill_between(f_range, atp_oxphos, atp_total, color='#4ade80', alpha=0.6, label='Glycolytic ATP')
ax.plot(f_range, atp_total, color='#e2e8f0', linewidth=2, label='Total ATP rate')
ax.axvline(x=f_glyc_default, color='#f87171', linestyle='--', linewidth=1.5, label=f'f = {f_glyc_default:.1f}')
ax.set_xlabel('Glycolytic fraction (f)', fontsize=11)
ax.set_ylabel('ATP production rate', fontsize=11)
ax.set_title('ATP Production vs Metabolic Strategy', fontsize=13)
ax.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1', fontsize=9)
ax.tick_params(colors='#cbd5e1')
ax.xaxis.label.set_color('#cbd5e1')
ax.yaxis.label.set_color('#cbd5e1')
ax.title.set_color('#e2e8f0')
for spine in ax.spines.values():
    spine.set_color('#334155')
ax.grid(True, alpha=0.15, color='#cbd5e1')

# Top-right: Lactate production
ax = axes[0, 1]
ax.set_facecolor('#1e293b')
ax.plot(f_range, lactate, color='#f87171', linewidth=2.5)
ax.fill_between(f_range, 0, lactate, color='#f87171', alpha=0.2)
ax.axvline(x=f_glyc_default, color='#fbbf24', linestyle='--', linewidth=1.5)
ax.set_xlabel('Glycolytic fraction (f)', fontsize=11)
ax.set_ylabel('Lactate production rate', fontsize=11)
ax.set_title('Lactate Production (Warburg Hallmark)', fontsize=13)
ax.tick_params(colors='#cbd5e1')
ax.xaxis.label.set_color('#cbd5e1')
ax.yaxis.label.set_color('#cbd5e1')
ax.title.set_color('#e2e8f0')
for spine in ax.spines.values():
    spine.set_color('#334155')
ax.grid(True, alpha=0.15, color='#cbd5e1')

# Bottom-left: Efficiency vs rate tradeoff
ax = axes[1, 0]
ax.set_facecolor('#1e293b')
efficiency = atp_total / (J_glc) if J_glc > 0 else f_range * 0
rate_advantage = atp_total / atp_total[0] * 100 if atp_total[0] > 0 else f_range * 0

ax.plot(f_range, efficiency, color='#c084fc', linewidth=2.5, label='ATP/glucose')
ax.set_xlabel('Glycolytic fraction (f)', fontsize=11)
ax.set_ylabel('ATP yield per glucose', fontsize=11)
ax.set_title('Efficiency: ATP per Glucose Consumed', fontsize=13)
ax.axvline(x=f_glyc_default, color='#f87171', linestyle='--', linewidth=1.5)
ax.tick_params(colors='#cbd5e1')
ax.xaxis.label.set_color('#cbd5e1')
ax.yaxis.label.set_color('#cbd5e1')
ax.title.set_color('#e2e8f0')
for spine in ax.spines.values():
    spine.set_color('#334155')
ax.grid(True, alpha=0.15, color='#cbd5e1')

# Bottom-right: Compare normal cell vs Warburg at different glucose uptakes
ax = axes[1, 1]
ax.set_facecolor('#1e293b')
glc_range = np.linspace(1, 30, 100)

# Normal cell: f=0.05, standard uptake
atp_normal = 2 * 0.05 * glc_range + 36 * 0.95 * glc_range * O2_availability / (K_O2 + O2_availability)

# Warburg (cancer): f=0.85, but can uptake much more glucose (upregulated GLUT1)
warburg_uptake = glc_range * 3  # 3x glucose uptake due to GLUT1 overexpression
atp_warburg = 2 * 0.85 * warburg_uptake + 36 * 0.15 * warburg_uptake * O2_availability / (K_O2 + O2_availability)

ax.plot(glc_range, atp_normal, color='#4ade80', linewidth=2.5, label='Normal cell (f=0.05, 1x uptake)')
ax.plot(glc_range, atp_warburg, color='#f87171', linewidth=2.5, label='Cancer cell (f=0.85, 3x uptake)')
ax.set_xlabel('Ambient glucose (relative)', fontsize=11)
ax.set_ylabel('ATP production rate', fontsize=11)
ax.set_title('Warburg Paradox: Rate vs Efficiency', fontsize=13)
ax.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1', fontsize=9)
ax.tick_params(colors='#cbd5e1')
ax.xaxis.label.set_color('#cbd5e1')
ax.yaxis.label.set_color('#cbd5e1')
ax.title.set_color('#e2e8f0')
for spine in ax.spines.values():
    spine.set_color('#334155')
ax.grid(True, alpha=0.15, color='#cbd5e1')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a', edgecolor='none')

# Compute values at default f
atp_at_f = 2 * f_glyc_default * J_glc + 36 * (1-f_glyc_default) * J_glc * O2_availability / (K_O2 + O2_availability)
eff_at_f = atp_at_f / J_glc if J_glc > 0 else 0
lac_at_f = 2 * f_glyc_default * J_glc

print("Warburg Effect Simulation")
print("=" * 45)
print(f"Glucose uptake: {J_glc:.1f} (a.u.)")
print(f"O2 availability: {O2_availability:.2f}")
print(f"Glycolytic fraction: {f_glyc_default:.2f}")
print()
print(f"At f = {f_glyc_default:.2f}:")
print(f"  Total ATP rate: {atp_at_f:.2f}")
print(f"  ATP/glucose: {eff_at_f:.2f}")
print(f"  Lactate rate: {lac_at_f:.2f}")
print()
print(f"At f = 0 (pure OxPhos):")
atp_0 = 36 * J_glc * O2_availability / (K_O2 + O2_availability)
print(f"  Total ATP rate: {atp_0:.2f}")
print(f"  ATP/glucose: {atp_0/J_glc:.2f}")
print()
print(f"At f = 1 (pure glycolysis):")
print(f"  Total ATP rate: {2*J_glc:.2f}")
print(f"  ATP/glucose: 2.00")
print(f"  Lactate rate: {2*J_glc:.2f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Understanding the Results

The top-left panel reveals the key insight: as the glycolytic fraction increases, total ATP production rate decreases because glycolysis is 18 times less efficient per glucose than oxidative phosphorylation. However, the bottom-right panel shows the resolution of the Warburg paradox: cancer cells compensate by massively upregulating glucose uptake (via GLUT1), so despite low efficiency, their absolute ATP production rate can exceed that of normal cells.

The high lactate production (top-right) is a clinical hallmark exploited in FDG-PET imaging, where the glucose analog fluorodeoxyglucose accumulates preferentially in tumors with high glycolytic rates. Try reducing O2 availability to see how hypoxia forces even normal cells toward glycolysis, mimicking the core of the tumor microenvironment.