Module 3: Microbial Biochemistry & Soil Ecosystems

Ecological Biochemistry & Biodiversity

1. Nitrogen Fixation

Biological nitrogen fixation is one of the most energetically demanding biochemical reactions in nature. The enzyme nitrogenase catalyzes the conversion of atmospheric dinitrogen (\(\text{N}_2\)) into bioavailable ammonia, a reaction that overcomes the formidable triple bond of N\(_2\) (bond energy\(\approx 945 \text{ kJ/mol}\)).

\[ \text{N}_2 + 8\text{H}^+ + 8e^- + 16\text{ATP} \longrightarrow 2\text{NH}_3 + \text{H}_2 + 16\text{ADP} + 16\text{P}_i \]

The overall nitrogenase reaction

Nitrogenase Mechanism: Two-Component System

Nitrogenase consists of two metalloprotein components that work in concert through a sophisticated electron transfer chain:

Fe-protein (dinitrogenase reductase) — a homodimer containing a single [4Fe-4S] cluster. It hydrolyzes ATP and transfers electrons one at a time to the MoFe-protein. Each electron transfer requires 2 ATP molecules:\[ \text{Fe-protein}_{\text{red}} + 2\text{ATP} \rightarrow \text{Fe-protein}_{\text{ox}} + 2\text{ADP} + 2\text{P}_i + e^- \]
MoFe-protein (dinitrogenase) — an \(\alpha_2\beta_2\) heterotetramer containing two unique metal clusters: the P-cluster ([8Fe-7S]) which receives electrons from Fe-protein, and the FeMo-cofactor (FeMoco, [7Fe-9S-Mo-C-homocitrate]) where N\(_2\) is actually reduced.

The electron transfer pathway:

\[ \text{Ferredoxin}_{\text{red}} \xrightarrow{} \text{Fe-protein [4Fe-4S]} \xrightarrow{2\text{ATP}} \text{P-cluster [8Fe-7S]} \xrightarrow{} \text{FeMoco} \xrightarrow{} \text{N}_2 \]

Energetic Cost Derivation

Why does nitrogen fixation require 16 ATP — making it one of the most expensive biological reactions? The answer involves both thermodynamic and kinetic factors.

The thermodynamic minimum (ignoring kinetic barriers):

\[ \Delta G°' = -33.5 \text{ kJ/mol (at pH 7, for N}_2 + 3\text{H}_2 \rightarrow 2\text{NH}_3) \]

This is actually exergonic! So why the ATP cost?

Kinetic barrier: The N\(\equiv\)N triple bond has an activation energy \(E_a \approx 945\) kJ/mol. Nitrogenase overcomes this by:

Binding N\(_2\) to FeMoco, weakening the triple bond through back-donation
Delivering electrons one at a time (8 total — 6 for reduction, 2 “wasted” on obligate H\(_2\) evolution)
Each electron transfer uses 2 ATP for conformational gating (ensuring unidirectional flow)

\[ \text{Total ATP} = 8 \text{ electrons} \times 2 \text{ ATP/electron} = 16 \text{ ATP} \]

\[ \text{Energy cost} = 16 \times 30.5 \text{ kJ/mol} = 488 \text{ kJ/mol N}_2 \]

Haber-Bosch Comparison

The industrial Haber-Bosch process (\(\text{N}_2 + 3\text{H}_2 \xrightarrow{400\text{-}500°C, 150\text{-}300\text{ atm}} 2\text{NH}_3\)) uses an iron catalyst at extreme conditions. While thermodynamically the same reaction, the industrial process consumes \(\sim 485\) kJ/mol N\(_2\) (including H\(_2\) production from methane steam reforming), comparable to the biological cost but achieved at ambient temperature and pressure by the enzyme. Haber-Bosch accounts for ~1-2% of global energy consumption and ~3-5% of natural gas production.

Rhizobium-Legume Symbiosis & Nod Factors

The most ecologically significant nitrogen fixation occurs in the mutualistic symbiosis between rhizobial bacteria and leguminous plants. This partnership involves an elaborate molecular dialogue:

Plant signal: Root exudates release flavonoids (e.g., luteolin, daidzein) that are specific to each legume species
Bacterial response: Flavonoids activate NodD transcription factors, inducing nod gene expression
Nod factors: Lipochitooligosaccharides (LCOs) — modified chitin oligomers with specific acyl chains — are secreted back to the plant
Root hair curling: Nod factors trigger Ca\(^{2+}\) spiking in root hair cells, causing them to curl and trap rhizobia in “infection threads”
Nodule formation: Cortical cell divisions create the root nodule — a specialized organ with leghemoglobin maintaining O\(_2\) at \(\sim 10\) nM (nitrogenase is irreversibly inhibited by O\(_2\))

The oxygen paradox — leghemoglobin equilibrium:

\[ \text{Lb} + \text{O}_2 \rightleftharpoons \text{LbO}_2 \quad K_d \approx 10 \text{ nM} \]

This maintains free O\(_2\) at ~10 nM while still delivering O\(_2\) to bacteroid respiration (needed for ATP production to fuel nitrogenase)

Global biological nitrogen fixation contributes approximately 140 Tg N/year — comparable to the 120 Tg N/year from Haber-Bosch. Together they have doubled the rate of nitrogen entering terrestrial ecosystems, with profound consequences for eutrophication, greenhouse gas emissions (N\(_2\)O), and biodiversity.

2. Decomposition Biochemistry

Decomposition of organic matter is mediated by extracellular enzymes secreted by soil microorganisms. Unlike intracellular enzymes that operate in controlled environments, extracellular enzymes face unique challenges: dilution in the soil matrix, adsorption to mineral surfaces, and substrate heterogeneity.

Key Decomposition Enzymes

Cellulase Complex

Endoglucanases cleave internal \(\beta\text{-1,4}\) glycosidic bonds; cellobiohydrolases attack chain ends; \(\beta\)-glucosidases hydrolyze cellobiose to glucose. Cellulose constitutes ~40-50% of plant biomass.

Lignin Peroxidase

White-rot fungi (e.g., Phanerochaete chrysosporium) produce LiP, MnP, and laccase. These generate radical species that non-specifically oxidize lignin's aromatic rings. Lignin is the second most abundant biopolymer (~25% of plant biomass).

Chitinase

Hydrolyzes chitin (\(\beta\text{-1,4}\)-linked N-acetylglucosamine), the main component of fungal cell walls and arthropod exoskeletons. Critical for nutrient recycling in soil and for fungal competition/antagonism.

Michaelis-Menten for Extracellular Enzymes

Extracellular enzymes follow modified Michaelis-Menten kinetics, but with important complications from the soil environment:

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

In soils, substrate limitation is common (\([S] \ll K_m\)), simplifying to first-order kinetics:

\[ v \approx \frac{V_{\max}}{K_m} \cdot [S] = k_{\text{cat/eff}} \cdot [S] \]

The apparent \(K_m\) in soil is modified by enzyme adsorption to mineral surfaces:

\[ K_m^{\text{app}} = K_m \cdot \left(1 + \frac{[E_{\text{bound}}]}{[E_{\text{free}}]}\right) = K_m \cdot \left(1 + K_{\text{ads}} \cdot [M]\right) \]

where \(K_{\text{ads}}\) is the adsorption coefficient and \([M]\) is mineral surface concentration

Decomposition Rate Model

The overall decomposition rate is modulated by environmental factors:

\[ k = k_{\max} \cdot f(T) \cdot f(\theta) \cdot f(Q) \]

Temperature function — Q\(_{10}\) model with optimum:

\[ f(T) = Q_{10}^{(T - T_{\text{ref}})/10} \cdot \exp\left(-\frac{(T - T_{\text{opt}})^2}{\sigma_T^2}\right) \]

Moisture function — Gaussian with optimum at ~50% saturation:

\[ f(\theta) = \exp\left(-\frac{(\theta - \theta_{\text{opt}})^2}{\sigma_\theta^2}\right) \]

Substrate quality — varies from labile (glucose, \(f = 1.0\)) to recalcitrant (lignin, \(f \approx 0.05\)):

\[ f(Q) = f_{\text{labile}} \cdot x_{\text{labile}} + f_{\text{cellulose}} \cdot x_{\text{cellulose}} + f_{\text{lignin}} \cdot x_{\text{lignin}} \]

Century Model: Three Carbon Pools

The CENTURY model (Parton et al., 1987) partitions soil organic matter into three pools with dramatically different turnover times:

Pool	Turnover Time	k (yr\(^{-1}\))	Composition
Active	1-5 years	0.5	Microbial biomass, labile metabolites
Slow	20-50 years	0.02	Partially decomposed plant material, resistant cell walls
Passive	200-1500 years	0.001	Humified material, mineral-associated organic matter

The coupled differential equations:

\[ \frac{dC_A}{dt} = I_A - k_A C_A \]

\[ \frac{dC_S}{dt} = I_S + f_{AS} k_A C_A - k_S C_S \]

\[ \frac{dC_P}{dt} = I_P + f_{AP} k_A C_A + f_{SP} k_S C_S - k_P C_P \]

where \(I\) = litter input, \(f_{ij}\) = transfer fraction from pool \(i\) to pool \(j\)

3. Soil Microbiome Metabolomics

A single gram of soil contains approximately \(10^9\) bacterial cells, representing\(10^3\) to \(10^4\) species — making soil the most microbiologically diverse habitat on Earth. Metagenomic sequencing reveals that the functional capacity of these communities is far more conserved than their taxonomic composition.

Functional Redundancy

Functional redundancy refers to the phenomenon where many taxonomically distinct species can perform the same metabolic function. For example, cellulose degradation is carried out by hundreds of bacterial and fungal species, each encoding similar cellulase gene families. This redundancy provides ecosystem resilience — loss of one species is compensated by others.

Quantifying redundancy with the redundancy ratio:

\[ R = \frac{S_{\text{species}}}{F_{\text{functions}}} \]

Typical soil values: \(R \approx 5\text{-}20\), meaning 5-20 species share each functional role. High redundancy \(\Rightarrow\) high resilience to species loss.

Shannon Diversity: Taxonomic vs Functional

Shannon diversity index quantifies both richness and evenness:

\[ H' = -\sum_{i=1}^{S} p_i \ln p_i \]

where \(p_i\) is the relative abundance of species/function \(i\)

Deriving the relationship between taxonomic and functional diversity:

If each functional group \(j\) contains \(n_j\) redundant species:

\[ H'_{\text{tax}} = H'_{\text{func}} + \sum_{j=1}^{F} q_j \cdot H'_{\text{within},j} \]

where \(q_j\) is the total abundance of functional group \(j\) and\(H'_{\text{within},j}\) is the within-group species diversity. This decomposition shows that:

\[ H'_{\text{tax}} \geq H'_{\text{func}} \]

Taxonomic diversity always exceeds functional diversity when redundancy exists. The gap\(\Delta H = H'_{\text{tax}} - H'_{\text{func}}\) quantifies the total functional redundancy in the community. In typical soils, \(H'_{\text{func}} / H'_{\text{tax}} \approx 0.4\text{-}0.7\).

Why Functional Diversity Matters More

Ecosystem process rates (decomposition, nutrient cycling, greenhouse gas flux) correlate more strongly with functional gene diversity than with species counts. Key examples:

Nitrogen cycling genes (nifH, amoA, nirS/nirK, nosZ) — the presence and expression of these genes, not the identity of the organisms carrying them, determines N\(_2\)O emissions
Methanogenesis genes (mcrA) — methane production depends on the functional capacity for methanogenesis, regardless of which archaea are present
Phosphorus solubilization — multiple bacterial lineages produce phosphatases and organic acids; the total enzymatic capacity determines plant P availability

This insight has transformed soil microbiology from a taxonomy-centered discipline to a function-centered one, with metatranscriptomics and metaproteomics revealing active metabolic pathways rather than mere species lists.

4. Carbon Cycling Enzymes

The global carbon cycle is driven by a handful of enzymes that collectively process ~120 Gt C/year through terrestrial photosynthesis and respiration. Understanding these enzymes is key to predicting carbon cycle feedbacks under climate change.

RuBisCO: The Most Abundant Enzyme on Earth

Ribulose-1,5-bisphosphate carboxylase/oxygenase (RuBisCO) catalyzes the first step of carbon fixation in the Calvin cycle. With an estimated 0.7 Gt on Earth, it constitutes ~50% of leaf protein and ~25% of leaf nitrogen.

Carboxylation reaction (desired):

\[ \text{RuBP} + \text{CO}_2 \xrightarrow{\text{RuBisCO}} 2 \times \text{3-PGA} \]

Oxygenation reaction (wasteful photorespiration):

\[ \text{RuBP} + \text{O}_2 \xrightarrow{\text{RuBisCO}} \text{3-PGA} + \text{2-PG (glycolate)} \]

The specificity factor determines CO\(_2\)/O\(_2\) selectivity:

\[ S_{\text{C/O}} = \frac{V_c K_O}{V_O K_C} \approx 80\text{-}100 \text{ (for C3 plants)} \]

RuBisCO's catalytic rate is remarkably slow (~3 s\(^{-1}\)), hence the need for such enormous quantities

Carbonic Anhydrase

One of the fastest enzymes known (\(k_{\text{cat}} \approx 10^6\) s\(^{-1}\)), carbonic anhydrase catalyzes CO\(_2\) hydration:

\[ \text{CO}_2 + \text{H}_2\text{O} \rightleftharpoons \text{H}_2\text{CO}_3 \rightleftharpoons \text{H}^+ + \text{HCO}_3^- \]

Essential for CO\(_2\) transport in blood, photosynthetic CO\(_2\) concentration, and ocean carbonate chemistry

Methane Monooxygenase (MMO)

Methanotrophs (CH\(_4\)-oxidizing bacteria) use MMO to convert methane to methanol:

\[ \text{CH}_4 + \text{O}_2 + \text{NAD(P)H} + \text{H}^+ \xrightarrow{\text{MMO}} \text{CH}_3\text{OH} + \text{NAD(P)}^+ + \text{H}_2\text{O} \]

Methanotrophs consume ~30 Tg CH\(_4\)/year, serving as a biological methane sink

Global Carbon Flux Budget Derivation

The carbon budget of an ecosystem is described by a hierarchy of production terms:

\[ \text{GPP} = \text{Total CO}_2 \text{ fixed by photosynthesis} \approx 120 \text{ Gt C/yr (terrestrial)} \]

\[ \text{NPP} = \text{GPP} - R_a \approx 60 \text{ Gt C/yr} \]

\[ \text{NEP} = \text{NPP} - R_h = \text{GPP} - R_a - R_h \approx 1\text{-}3 \text{ Gt C/yr} \]

\[ \text{NBP} = \text{NEP} - D_{\text{fire}} - D_{\text{harvest}} - D_{\text{erosion}} \approx 0.5\text{-}2 \text{ Gt C/yr} \]

Where:

\(R_a\) = autotrophic respiration (plant mitochondrial respiration) \(\approx 60\) Gt C/yr
\(R_h\) = heterotrophic respiration (decomposers, soil microbes) \(\approx 57\text{-}59\) Gt C/yr
\(D\) = disturbance losses (fire, harvest, erosion)
NBP = Net Biome Production — the actual carbon sequestration rate

The terrestrial biosphere is currently a weak net carbon sink (NBP \(> 0\)), absorbing ~25% of anthropogenic CO\(_2\) emissions. However, warming-driven increases in \(R_h\)(decomposition) may eventually exceed photosynthetic gains, potentially converting terrestrial ecosystems from sinks to sources — a critical tipping point in climate science.

5. Soil Carbon Cycle Diagram

The following diagram illustrates the flow of carbon through soil ecosystems, from atmospheric CO\(_2\) fixation through plant tissues to microbial decomposition pools and back:

6. Computational Simulations

The following simulations model: (1) decomposition rate as a function of temperature, moisture, and substrate quality; (2) nitrogen fixation energetics comparing biological and industrial processes; (3) Century-like 3-pool soil carbon dynamics; and (4) functional vs taxonomic diversity showing redundancy patterns.

Python

script.py235 lines

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

# ══════════════════════════════════════════════════════════════════════
# Panel 1: Decomposition rate as function of Temperature x Moisture x Quality
# ══════════════════════════════════════════════════════════════════════
T = np.linspace(0, 50, 200)         # temperature (°C)
moisture = np.linspace(0, 1, 200)    # volumetric water content (fraction)

# Temperature response (Q10 with optimum)
def f_temp(T, Topt=25, Q10=2.0):
    """Bell-shaped temperature response"""
    k_ref = 1.0
    rising = k_ref * Q10 ** ((T - 20) / 10)
    decline = np.exp(-((T - Topt) / 12) ** 2)
    return rising * decline / (k_ref * Q10 ** ((Topt - 20) / 10))

# Moisture response (optimum around 0.4-0.6)
def f_moisture(w, w_opt=0.5, sigma=0.2):
    return np.exp(-((w - w_opt) / sigma) ** 2)

# Substrate quality response
qualities = {'Glucose (labile)': 1.0, 'Cellulose': 0.4, 'Lignin': 0.05, 'Humus': 0.01}

f_T = f_temp(T)
f_M = f_moisture(moisture)

# ══════════════════════════════════════════════════════════════════════
# Panel 2: Nitrogen Fixation Energetics
# ══════════════════════════════════════════════════════════════════════
# ATP cost analysis: N2 + 8H+ + 8e- + 16ATP -> 2NH3 + H2 + 16ADP
# Haber-Bosch: N2 + 3H2 -> 2NH3 at 400-500°C, 150-300 atm
atp_per_n2 = 16
energy_per_atp = 30.5  # kJ/mol (standard)
total_bio_energy = atp_per_n2 * energy_per_atp  # kJ per mol N2
haber_energy = 485  # kJ per mol N2 (industrial, including H2 production)

# Compare with other N sources
n_sources = ['Biological\nfixation', 'Haber-Bosch\n(total)', 'Lightning\nfixation', 'Ammonification\n(recycling)']
energies = [total_bio_energy, haber_energy, 2500, 15]
global_flux = [140, 120, 5, 400]  # Tg N/year

# ══════════════════════════════════════════════════════════════════════
# Panel 3: Century Model — 3-pool soil carbon simulation
# ══════════════════════════════════════════════════════════════════════
dt = 0.01   # years
t_max = 200
t = np.arange(0, t_max, dt)
n = len(t)

C_active = np.zeros(n)     # active pool (turnover ~1-5 yr)
C_slow = np.zeros(n)       # slow pool (turnover ~20-50 yr)
C_passive = np.zeros(n)    # passive pool (turnover ~200-1500 yr)
C_total = np.zeros(n)

# Initial conditions (g C / m^2)
C_active[0] = 200
C_slow[0] = 3000
C_passive[0] = 8000

# Rate constants (yr^-1)
k_active = 0.5      # fast turnover
k_slow = 0.02       # moderate turnover
k_passive = 0.001   # very slow turnover

# Transfer fractions
f_as = 0.4    # active -> slow
f_ap = 0.05   # active -> passive
f_sp = 0.03   # slow -> passive

# Litter input (g C / m^2 / yr)
litter_input = 400

# Partitioning of litter
f_litter_active = 0.55
f_litter_slow = 0.40
f_litter_passive = 0.05

for i in range(n - 1):
    decay_a = k_active * C_active[i]
    decay_s = k_slow * C_slow[i]
    decay_p = k_passive * C_passive[i]

C_active[i+1] = C_active[i] + dt * (
        litter_input * f_litter_active - decay_a + 0)
    C_slow[i+1] = C_slow[i] + dt * (
        litter_input * f_litter_slow + f_as * decay_a - decay_s)
    C_passive[i+1] = C_passive[i] + dt * (
        litter_input * f_litter_passive + f_ap * decay_a + f_sp * decay_s - decay_p)

C_total = C_active + C_slow + C_passive

# ══════════════════════════════════════════════════════════════════════
# Panel 4: Shannon Diversity — Functional vs Taxonomic
# ══════════════════════════════════════════════════════════════════════
np.random.seed(42)
n_communities = 50
taxonomic_H = np.zeros(n_communities)
functional_H = np.zeros(n_communities)

for c in range(n_communities):
    n_species = np.random.randint(10, 200)
    # Species abundances (log-normal)
    abundances = np.random.lognormal(0, 1.5, n_species)
    p_tax = abundances / abundances.sum()
    taxonomic_H[c] = -np.sum(p_tax * np.log(p_tax))

# Functional genes — fewer unique functions due to redundancy
    n_functions = max(5, int(n_species * np.random.uniform(0.1, 0.5)))
    func_abundances = np.random.lognormal(0, 1.0, n_functions)
    p_func = func_abundances / func_abundances.sum()
    functional_H[c] = -np.sum(p_func * np.log(p_func))

# ══════════════════════════════════════════════════════════════════════
# Plotting
# ══════════════════════════════════════════════════════════════════════
fig, axes = plt.subplots(2, 2, figsize=(15, 12))
fig.patch.set_facecolor('#0a0a1a')

# Plot 1: Decomposition rate components
ax1 = axes[0, 0]
ax1.plot(T, f_T, color='#2dd4bf', linewidth=2.5, label='f(Temperature)')
ax1_twin = ax1.twinx()
ax1_twin.plot(moisture * 100, f_M, color='#60a5fa', linewidth=2.5, label='f(Moisture)')
ax1.set_xlabel('Temperature (°C) / Moisture (%)', color='white', fontsize=11)
ax1.set_ylabel('Temperature response', color='#2dd4bf', fontsize=11)
ax1_twin.set_ylabel('Moisture response', color='#60a5fa', fontsize=11)
ax1.set_title('Decomposition Rate Modifiers\nk = k_max · f(T) · f(moisture) · f(quality)', color='white', fontsize=12, fontweight='bold')

# Add substrate quality bars inset
inset_ax = fig.add_axes([0.15, 0.72, 0.12, 0.12])
inset_ax.barh(list(qualities.keys()), list(qualities.values()), color='#34d399', alpha=0.8)
inset_ax.set_xlabel('f(quality)', color='white', fontsize=8)
inset_ax.set_facecolor('#0a0a1a')
inset_ax.tick_params(colors='white', labelsize=7)
for sp in inset_ax.spines.values(): sp.set_color('#2dd4bf')

ax1.set_facecolor('#0a0a1a')
ax1.tick_params(colors='white')
ax1_twin.tick_params(colors='white')
ax1.grid(True, alpha=0.15, color='#2dd4bf')
for sp in ax1.spines.values(): sp.set_color('#2dd4bf')
for sp in ax1_twin.spines.values(): sp.set_color('#2dd4bf')

lines1, labels1 = ax1.get_legend_handles_labels()
lines2, labels2 = ax1_twin.get_legend_handles_labels()
ax1.legend(lines1 + lines2, labels1 + labels2, fontsize=9, facecolor='#1a1a2e', edgecolor='#2dd4bf', labelcolor='white')

# Plot 2: Nitrogen fixation energetics
ax2 = axes[0, 1]
colors = ['#2dd4bf', '#f87171', '#fbbf24', '#a78bfa']
bars = ax2.bar(n_sources, energies, color=colors, alpha=0.85, edgecolor='white', linewidth=0.8)
for bar, e in zip(bars, energies):
    ax2.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 20,
             f'{e:.0f}', ha='center', va='bottom', color='white', fontsize=10, fontweight='bold')
ax2.set_ylabel('Energy cost (kJ / mol N₂)', color='white', fontsize=11)
ax2.set_title('Nitrogen Fixation Energy Cost Comparison', color='white', fontsize=12, fontweight='bold')
ax2.set_facecolor('#0a0a1a')
ax2.tick_params(colors='white')
ax2.grid(True, alpha=0.15, color='#2dd4bf', axis='y')
for sp in ax2.spines.values(): sp.set_color('#2dd4bf')

# Add global flux as secondary axis
ax2b = ax2.twinx()
ax2b.plot(n_sources, global_flux, 'D-', color='#fbbf24', markersize=8, linewidth=2, label='Global flux (Tg N/yr)')
ax2b.set_ylabel('Global N flux (Tg N/yr)', color='#fbbf24', fontsize=11)
ax2b.tick_params(colors='white')
ax2b.legend(fontsize=9, facecolor='#1a1a2e', edgecolor='#fbbf24', labelcolor='white', loc='upper right')
for sp in ax2b.spines.values(): sp.set_color('#2dd4bf')

# Plot 3: Century model carbon pools
ax3 = axes[1, 0]
ax3.fill_between(t, 0, C_active, alpha=0.6, color='#2dd4bf', label=f'Active (k={k_active}/yr)')
ax3.fill_between(t, C_active, C_active + C_slow, alpha=0.6, color='#60a5fa', label=f'Slow (k={k_slow}/yr)')
ax3.fill_between(t, C_active + C_slow, C_total, alpha=0.6, color='#a78bfa', label=f'Passive (k={k_passive}/yr)')
ax3.plot(t, C_total, color='white', linewidth=2, label='Total SOC')
ax3.set_xlabel('Time (years)', color='white', fontsize=11)
ax3.set_ylabel('Soil Organic Carbon (g C/m²)', color='white', fontsize=11)
ax3.set_title('Century Model: Soil Carbon Pool Dynamics', color='white', fontsize=12, fontweight='bold')
ax3.legend(fontsize=9, facecolor='#1a1a2e', edgecolor='#2dd4bf', labelcolor='white')
ax3.set_facecolor('#0a0a1a')
ax3.tick_params(colors='white')
ax3.grid(True, alpha=0.15, color='#2dd4bf')
for sp in ax3.spines.values(): sp.set_color('#2dd4bf')

# Plot 4: Functional vs Taxonomic diversity
ax4 = axes[1, 1]
ax4.scatter(taxonomic_H, functional_H, c='#2dd4bf', alpha=0.7, s=60, edgecolors='white', linewidths=0.5)

# Fit line
z = np.polyfit(taxonomic_H, functional_H, 1)
p = np.poly1d(z)
x_fit = np.linspace(taxonomic_H.min(), taxonomic_H.max(), 100)
ax4.plot(x_fit, p(x_fit), '--', color='#f87171', linewidth=2,
         label=f'Linear fit (slope={z[0]:.2f})')
ax4.plot([0, 6], [0, 6], ':', color='#fbbf24', alpha=0.5, label='1:1 line')
ax4.set_xlabel('Taxonomic Shannon Diversity (H)', color='white', fontsize=11)
ax4.set_ylabel('Functional Shannon Diversity (H)', color='white', fontsize=11)
ax4.set_title('Functional vs Taxonomic Diversity\n(Functional Redundancy)', color='white', fontsize=12, fontweight='bold')
ax4.legend(fontsize=9, facecolor='#1a1a2e', edgecolor='#2dd4bf', labelcolor='white')
ax4.set_facecolor('#0a0a1a')
ax4.tick_params(colors='white')
ax4.grid(True, alpha=0.15, color='#2dd4bf')
for sp in ax4.spines.values(): sp.set_color('#2dd4bf')

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

# ── Numerical Output ──
print("=== Microbial Soil Biochemistry Simulations ===\n")

print("--- Decomposition Rate Analysis ---")
print(f"Peak temperature response at Topt = 25°C")
print(f"Q10 factor = 2.0 (rate doubles per 10°C rise)")
print(f"Optimal moisture ~50% volumetric water content\n")

print("--- Nitrogen Fixation Energetics ---")
print(f"Biological: {atp_per_n2} ATP x {energy_per_atp} kJ/mol = {total_bio_energy} kJ/mol N2")
print(f"Haber-Bosch: ~{haber_energy} kJ/mol N2 (industrial)")
print(f"Energy ratio (Haber/Bio): {haber_energy/total_bio_energy:.1f}x\n")

print("--- Century Model Steady State ---")
print(f"Active pool: {C_active[-1]:.0f} g C/m²")
print(f"Slow pool: {C_slow[-1]:.0f} g C/m²")
print(f"Passive pool: {C_passive[-1]:.0f} g C/m²")
print(f"Total SOC: {C_total[-1]:.0f} g C/m²\n")

print("--- Diversity Analysis ---")
print(f"Mean taxonomic H: {taxonomic_H.mean():.2f}")
print(f"Mean functional H: {functional_H.mean():.2f}")
print(f"Functional/Taxonomic ratio: {(functional_H/taxonomic_H).mean():.2f}")
print(f"This demonstrates functional redundancy: many species share the same metabolic roles.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

7. Denitrification & the Nitrogen Cascade

Denitrification is the stepwise microbial reduction of nitrate to dinitrogen gas, returning reactive nitrogen to the atmosphere. Each step is catalyzed by a distinct metalloenzyme, and the pathway represents a major source of the potent greenhouse gas N\(_2\)O.

\[ \text{NO}_3^- \xrightarrow{\text{Nar/Nap}} \text{NO}_2^- \xrightarrow{\text{NirS/NirK}} \text{NO} \xrightarrow{\text{Nor}} \text{N}_2\text{O} \xrightarrow{\text{NosZ}} \text{N}_2 \]

The complete denitrification pathway

Metalloenzymes of the Nitrogen Cascade

Nitrate Reductase (Nar/Nap)

Contains a molybdenum (Mo) cofactor at the active site (Mo-bis-molybdopterin guanine dinucleotide). Nar is membrane-bound and coupled to the proton motive force; Nap is periplasmic. Catalyzes:\(\text{NO}_3^- + 2\text{H}^+ + 2e^- \rightarrow \text{NO}_2^- + \text{H}_2\text{O}\)

Nitrite Reductase (NirS/NirK)

Two non-homologous forms: NirS contains heme cd\(_1\) (Fe); NirK contains copper (Cu) at the type 1 and type 2 sites. This is the committed step — once NO is produced, the cell is committed to denitrification:\(\text{NO}_2^- + 2\text{H}^+ + e^- \rightarrow \text{NO} + \text{H}_2\text{O}\)

Nitric Oxide Reductase (Nor)

Contains a binuclear iron (Fe) center (heme b\(_3\)-Fe\(_B\)). Related to cytochrome c oxidase. Detoxifies the cytotoxic NO radical:\(2\text{NO} + 2\text{H}^+ + 2e^- \rightarrow \text{N}_2\text{O} + \text{H}_2\text{O}\)

Nitrous Oxide Reductase (NosZ)

Contains the unique Cu\(_Z\) cluster — a tetranuclear copper-sulfide center [4Cu:2S]. The only known enzyme that reduces N\(_2\)O. Many denitrifiers lack nosZ, making them net N\(_2\)O sources:\(\text{N}_2\text{O} + 2\text{H}^+ + 2e^- \rightarrow \text{N}_2 + \text{H}_2\text{O}\)

Gibbs Free Energy for Each Step

Each denitrification step is exergonic under standard conditions (pH 7), providing energy for anaerobic respiration:

\[ \text{NO}_3^- \rightarrow \text{NO}_2^-: \quad \Delta G°' = -163 \text{ kJ/mol (per 2e}^-) \]

\[ \text{NO}_2^- \rightarrow \text{NO}: \quad \Delta G°' = -73 \text{ kJ/mol (per e}^-) \]

\[ 2\text{NO} \rightarrow \text{N}_2\text{O}: \quad \Delta G°' = -306 \text{ kJ/mol (per 2e}^-) \]

\[ \text{N}_2\text{O} \rightarrow \text{N}_2: \quad \Delta G°' = -340 \text{ kJ/mol (per 2e}^-) \]

Total: \(\Delta G°' = -882\) kJ/mol for complete denitrification (\(\text{NO}_3^- \rightarrow \frac{1}{2}\text{N}_2\))

N\(_2\)O: A Potent Greenhouse Gas

Nitrous oxide has a Global Warming Potential (GWP) of 298 over a 100-year horizon — meaning 1 kg of N\(_2\)O traps as much heat as 298 kg of CO\(_2\). It is also the dominant ozone-depleting substance emitted today. Agricultural soils are the largest anthropogenic source (~4.2 Tg N\(_2\)O-N/yr), driven by nitrogen fertilizer application.

The Leaky Pipe Model

The “leaky pipe” model (Firestone & Davidson, 1989) treats the nitrogen cycle as a pipe where gaseous intermediates (NO, N\(_2\)O) leak out at each step. The fraction lost depends critically on oxygen availability:

The N\(_2\)O/(N\(_2\)O + N\(_2\)) ratio as a function of soil moisture (WFPS = water-filled pore space):

\[ r = \frac{\text{N}_2\text{O}}{\text{N}_2\text{O} + \text{N}_2} = \frac{1}{1 + \exp\!\left(k \cdot (\text{WFPS} - \text{WFPS}_{\text{crit}})\right)} \]

where:

At low WFPS (<60%): aerobic conditions dominate, nitrification is the main N\(_2\)O source, \(r \approx 0.9\text{-}1.0\)
At intermediate WFPS (60-80%): incomplete denitrification, high N\(_2\)O emissions, \(r \approx 0.3\text{-}0.7\)
At high WFPS (>80%): complete denitrification, N\(_2\) dominates, \(r \approx 0.0\text{-}0.2\)

Deriving O\(_2\) control of NosZ activity:

The Cu\(_Z\) center in NosZ is extremely oxygen-sensitive. The enzyme's activity depends on O\(_2\) concentration:

\[ v_{\text{NosZ}} = V_{\max} \cdot \frac{[\text{N}_2\text{O}]}{K_m + [\text{N}_2\text{O}]} \cdot \frac{K_i}{K_i + [\text{O}_2]} \]

where \(K_i \approx 0.3\) \(\mu\)M — NosZ is inhibited at O\(_2\) concentrations well below atmospheric. This is why partially waterlogged soils (intermediate O\(_2\)) produce the most N\(_2\)O: enough anoxia for denitrification to proceed, but too much O\(_2\) for the final step to N\(_2\).

8. Phosphorus Cycling & Limitation

Unlike nitrogen and carbon, phosphorus has no significant gaseous phase — it cycles entirely through the lithosphere, hydrosphere, and biosphere. This makes phosphorus fundamentally different from other macronutrients: the only new input to ecosystems comes from rock weathering, an extremely slow geological process.

Phosphatases: Releasing P from Organic Matter

In most soils, 30-80% of total P is in organic form, bound in phospholipids, nucleic acids, and phytic acid (inositol hexakisphosphate). Phosphatases are extracellular enzymes that hydrolyze ester bonds to release inorganic phosphate (P\(_i\)):

\[ \text{R-O-PO}_3^{2-} + \text{H}_2\text{O} \xrightarrow{\text{phosphatase}} \text{R-OH} + \text{HPO}_4^{2-} \]

Acid phosphatases (optimal pH 4-6): dominant in acidic soils, produced by plant roots and fungi. Alkaline phosphatases (optimal pH 8-10): produced by bacteria, dominant in calcareous soils. Expression is strongly upregulated under P limitation — a classic example of nutrient-responsive gene regulation.

The Redfield Ratio

Alfred Redfield (1934) discovered a remarkably constant elemental ratio in marine phytoplankton and dissolved nutrients:

\[ \text{C}:\text{N}:\text{P} = 106:16:1 \]

The Redfield ratio — a stoichiometric fingerprint of ocean biochemistry

This ratio arises from the average molecular composition of phytoplankton biomass:

\[ (\text{CH}_2\text{O})_{106}(\text{NH}_3)_{16}(\text{H}_3\text{PO}_4) + 138\text{O}_2 \rightarrow 106\text{CO}_2 + 122\text{H}_2\text{O} + 16\text{HNO}_3 + \text{H}_3\text{PO}_4 \]

The Redfield ratio has been remarkably stable across ocean basins and geological time, suggesting that phytoplankton community composition adjusts to reflect the nutrient supply ratio — a phenomenon called Redfield regulation.

Liebig's Law of the Minimum

Growth is limited by the nutrient in shortest supply relative to demand. Mathematically, this is expressed as:

\[ \mu = \mu_{\max} \cdot \min_i\!\left(\frac{S_i}{K_i + S_i}\right) \quad \text{for } i = \text{N, P, Fe, ...} \]

Growth rate \(\mu\) is determined by the most limiting nutrient

Each nutrient follows Monod (Michaelis-Menten) kinetics independently:

\[ f_N = \frac{[\text{N}]}{K_N + [\text{N}]}, \quad f_P = \frac{[\text{P}]}{K_P + [\text{P}]}, \quad f_{\text{Fe}} = \frac{[\text{Fe}]}{K_{\text{Fe}} + [\text{Fe}]} \]

Then: \(\mu = \mu_{\max} \cdot \min(f_N, f_P, f_{\text{Fe}})\)

Typical half-saturation constants: \(K_N \approx 0.5\text{-}2\) \(\mu\)M,\(K_P \approx 0.01\text{-}0.1\) \(\mu\)M,\(K_{\text{Fe}} \approx 0.01\text{-}0.1\) nM.

Why P is the Ultimate Limiting Nutrient in Freshwater

While nitrogen limits productivity in much of the ocean, phosphorus is the ultimate limiting nutrient in freshwater systems. The reasoning:

N deficiency can be overcome by N\(_2\)-fixing cyanobacteria — the atmosphere provides an inexhaustible N reservoir
P has no atmospheric reservoir and no biological fixation pathway — the only new P comes from rock weathering at ~0.1 Tg P/yr
On geological timescales, N-fixers expand until P becomes limiting, making P the ultimate constraint
Anthropogenic P loading (fertilizer, detergents) causes eutrophication: excess P \(\rightarrow\) algal blooms \(\rightarrow\) hypoxia \(\rightarrow\) dead zones

Global P budget:

Rock weathering input: ~10-20 Tg P/yr
Fertilizer application: ~20 Tg P/yr (mined phosphate rock — a non-renewable resource)
River transport to ocean: ~4 Tg P/yr (dissolved), ~20 Tg P/yr (particulate)
Ocean residence time: ~20,000-80,000 years
Peak phosphorus: global reserves may be depleted in 50-100 years at current extraction rates

9. Methanogens & the Methane Cycle

Methanogenesis is the terminal step in anaerobic decomposition, carried out exclusively by methanogenic archaea — obligate anaerobes that occupy the most reduced ecological niche on Earth. They produce approximately 1 Gt CH\(_4\)/yr globally, making biological methanogenesis the largest natural source of atmospheric methane.

Two Major Methanogenic Pathways

Aceticlastic Methanogenesis

Accounts for ~70% of biogenic methane. Acetate is cleaved by Methanosarcinaand Methanosaeta:

\[ \text{CH}_3\text{COOH} \rightarrow \text{CH}_4 + \text{CO}_2 \]

Hydrogenotrophic Methanogenesis

Accounts for ~30% of biogenic methane. Uses H\(_2\) as electron donor to reduce CO\(_2\) (most methanogen genera):

\[ \text{CO}_2 + 4\text{H}_2 \rightarrow \text{CH}_4 + 2\text{H}_2\text{O} \]

Thermodynamics of Methanogenic Pathways

Standard free energies under biological conditions (pH 7, 25\(°\)C):

\[ \text{Aceticlastic: } \text{CH}_3\text{COO}^- + \text{H}_2\text{O} \rightarrow \text{CH}_4 + \text{HCO}_3^- \quad \Delta G°' = -31 \text{ kJ/mol} \]

\[ \text{Hydrogenotrophic: } 4\text{H}_2 + \text{HCO}_3^- + \text{H}^+ \rightarrow \text{CH}_4 + 3\text{H}_2\text{O} \quad \Delta G°' = -136 \text{ kJ/mol} \]

Note the remarkably small \(\Delta G°'\) for the aceticlastic pathway — methanogens operate at the thermodynamic limit of life. Under actual conditions with low H\(_2\) partial pressures:

\[ \Delta G = \Delta G°' + RT \ln\!\left(\frac{p_{\text{CH}_4}}{p_{\text{H}_2}^4 \cdot p_{\text{CO}_2}}\right) \]

Hydrogenotrophic methanogens maintain H\(_2\) at \(\sim 10\) Pa — below this threshold, the reaction becomes endergonic and methanogenesis ceases. This creates the concept of a thermodynamic threshold for microbial metabolism.

Methanotrophs: The Biological Methane Sink

Methanotrophic bacteria oxidize CH\(_4\) using methane monooxygenase (MMO), consuming approximately 30 Tg CH\(_4\)/yr — acting as a critical biological filter that prevents much of the methane produced in soils and sediments from reaching the atmosphere.

\[ \text{CH}_4 + \text{O}_2 + \text{NAD(P)H} + \text{H}^+ \xrightarrow{\text{MMO}} \text{CH}_3\text{OH} + \text{NAD(P)}^+ + \text{H}_2\text{O} \]

Two forms exist: soluble MMO (sMMO, contains a di-iron center) and particulate MMO (pMMO, contains copper). pMMO is the dominant form in nature and is expressed when copper is available.

Wetland Methane Emissions Model

Wetlands are the largest natural methane source (~150-200 Tg CH\(_4\)/yr). Net emissions represent the balance between methanogenesis and methanotrophy:

\[ F_{\text{CH}_4} = P_{\text{CH}_4} \cdot (1 - f_{\text{ox}}) \cdot f_{\text{transport}} \]

where:

\(P_{\text{CH}_4}\) = gross methanogenesis rate (depends on T, substrate availability, redox)
\(f_{\text{ox}}\) = fraction oxidized by methanotrophs in the oxic zone (~60-90%)
\(f_{\text{transport}}\) = transport factor (diffusion, ebullition, plant-mediated transport via aerenchyma)

Temperature sensitivity of wetland CH\(_4\) emissions:

\[ P_{\text{CH}_4}(T) = P_{\text{ref}} \cdot Q_{10}^{(T - T_{\text{ref}})/10} \]

With \(Q_{10} \approx 3\text{-}5\) for methanogenesis (higher than most biogeochemical processes), wetland methane emissions are highly sensitive to warming — a potential positive climate feedback. Arctic permafrost thaw could release 50-100 Pg C as CH\(_4\)and CO\(_2\) over the coming century.

10. Advanced Simulations: N Cascade, P Limitation & Methane Balance

The following simulations model: (1) the complete denitrification nitrogen cascade with N\(_2\)O emissions as a function of soil moisture; (2) phosphorus limitation using the Droop (cell quota) model showing growth rate dependence on internal P stores; (3) methane production/consumption balance in wetlands across temperature.

Python

script.py222 lines

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

# ══════════════════════════════════════════════════════════════════════
# Panel 1: Nitrogen Cascade — N₂O Emissions vs Soil Moisture (WFPS)
# ══════════════════════════════════════════════════════════════════════
WFPS = np.linspace(0, 1, 300)  # Water-Filled Pore Space (fraction)

# Denitrification rate increases with WFPS (anaerobic)
denitrif_rate = np.where(WFPS < 0.5, 0.01, (WFPS - 0.5) / 0.5) ** 1.5

# Nitrification rate peaks at moderate WFPS
nitrif_rate = np.exp(-((WFPS - 0.55) / 0.15) ** 2)

# N₂O / (N₂O + N₂) ratio — decreases as WFPS increases (more complete denitrification)
r_N2O = 1.0 / (1.0 + np.exp(12 * (WFPS - 0.78)))

# Total N₂O emissions = nitrification N₂O + denitrification N₂O
N2O_nitrif = nitrif_rate * 0.01  # ~1% of nitrification lost as N₂O
N2O_denitrif = denitrif_rate * r_N2O * 0.8
N2O_total = N2O_nitrif + N2O_denitrif
N2_emission = denitrif_rate * (1 - r_N2O) * 0.8

# ══════════════════════════════════════════════════════════════════════
# Panel 2: Phosphorus Limitation — Droop Model (Cell Quota)
# ══════════════════════════════════════════════════════════════════════
dt = 0.01
t_max = 30
t = np.arange(0, t_max, dt)
n = len(t)

# State variables
B = np.zeros(n)    # Biomass (cells/mL)
S_P = np.zeros(n)  # External P (umol/L)
Q = np.zeros(n)    # Internal P quota (umol P / cell)

B[0] = 1e4
S_P[0] = 5.0    # initial external P
Q[0] = 2e-6     # initial cell quota

mu_max = 1.5     # max growth rate (day⁻¹)
Q_min = 5e-7     # minimum cell quota (subsistence)
V_max = 8e-6     # max P uptake rate (umol P / cell / day)
K_P = 0.5        # half-saturation for P uptake (umol/L)

for i in range(n - 1):
    # Droop growth: mu = mu_max * (1 - Q_min/Q)
    mu = mu_max * max(0, (1 - Q_min / max(Q[i], Q_min * 1.001)))
    # Monod uptake
    V = V_max * S_P[i] / (K_P + S_P[i])

B[i+1] = B[i] + dt * mu * B[i]
    S_P[i+1] = max(0, S_P[i] - dt * V * B[i])
    Q[i+1] = Q[i] + dt * (V - mu * Q[i])

if B[i+1] > 1e9: B[i+1] = 1e9  # cap

# ══════════════════════════════════════════════════════════════════════
# Panel 3: Methane Production/Consumption Balance in Wetlands
# ══════════════════════════════════════════════════════════════════════
T_wetland = np.linspace(-5, 35, 200)  # Temperature (°C)

# Methanogenesis rate (Q10 ~ 4)
P_ref = 50  # mg CH4 / m² / day at 15°C
Q10_meth = 4.0
P_CH4 = P_ref * Q10_meth ** ((T_wetland - 15) / 10)
P_CH4 = np.where(T_wetland < 0, P_CH4 * 0.05, P_CH4)  # frozen soil minimal

# Methanotrophy (oxidation fraction, ~60-90%)
f_ox_base = 0.7  # base oxidation fraction
Q10_ox = 2.0
methanotrophy = P_CH4 * f_ox_base * Q10_ox ** ((T_wetland - 15) / 10) / Q10_meth ** ((T_wetland - 15) / 10)
methanotrophy = np.minimum(methanotrophy, P_CH4 * 0.9)

# Net emission
net_CH4 = P_CH4 - methanotrophy

# Also compute for different water table depths
wt_shallow = net_CH4 * 1.3   # shallow water table = more emission
wt_deep = net_CH4 * 0.4      # deep water table = more oxidation

# ══════════════════════════════════════════════════════════════════════
# Panel 4: Gibbs Free Energy Cascade for Denitrification Steps
# ══════════════════════════════════════════════════════════════════════
steps = ['NO₃⁻\n→ NO₂⁻', 'NO₂⁻\n→ NO', '2NO\n→ N₂O', 'N₂O\n→ N₂']
dG_values = [-163, -73, -306, -340]
metals = ['Mo\n(Nar)', 'Cu/Fe\n(NirS/K)', 'Fe\n(Nor)', 'Cu₄S₂\n(NosZ)']
cumulative = np.cumsum([0] + dG_values)

# Plot 1: N₂O emissions vs WFPS
ax1 = axes[0, 0]
ax1.fill_between(WFPS * 100, 0, N2O_nitrif, alpha=0.4, color='#fbbf24', label='N₂O from nitrification')
ax1.fill_between(WFPS * 100, N2O_nitrif, N2O_total, alpha=0.4, color='#f87171', label='N₂O from denitrification')
ax1.plot(WFPS * 100, N2O_total, color='#f87171', linewidth=2.5, label='Total N₂O')
ax1.plot(WFPS * 100, N2_emission, color='#2dd4bf', linewidth=2, linestyle='--', label='N₂ emission')

ax1_twin = ax1.twinx()
ax1_twin.plot(WFPS * 100, r_N2O, color='#a78bfa', linewidth=2, linestyle=':', label='N₂O/(N₂O+N₂) ratio')
ax1_twin.set_ylabel('N₂O/(N₂O+N₂) ratio', color='#a78bfa', fontsize=11)
ax1_twin.tick_params(colors='white')
ax1_twin.set_ylim(-0.05, 1.1)
for sp in ax1_twin.spines.values(): sp.set_color('#2dd4bf')

ax1.set_xlabel('Water-Filled Pore Space (%)', color='white', fontsize=11)
ax1.set_ylabel('Relative Emission Rate', color='white', fontsize=11)
ax1.set_title('Leaky Pipe: N₂O Emissions vs Soil Moisture\n(GWP of N₂O = 298)', color='white', fontsize=12, fontweight='bold')
lines1, labels1 = ax1.get_legend_handles_labels()
lines2, labels2 = ax1_twin.get_legend_handles_labels()
ax1.legend(lines1 + lines2, labels1 + labels2, fontsize=8, facecolor='#1a1a2e', edgecolor='#2dd4bf', labelcolor='white', loc='upper left')
ax1.set_facecolor('#0a0a1a')
ax1.tick_params(colors='white')
ax1.grid(True, alpha=0.15, color='#2dd4bf')
for sp in ax1.spines.values(): sp.set_color('#2dd4bf')

# Plot 2: Droop model P limitation
ax2 = axes[0, 1]
color_B = '#2dd4bf'
color_P = '#f87171'
color_Q = '#fbbf24'

ln1 = ax2.plot(t, B / 1e4, color=color_B, linewidth=2.5, label='Biomass (×10⁴ cells/mL)')
ax2.set_ylabel('Biomass (×10⁴ cells/mL)', color=color_B, fontsize=11)

ax2b = ax2.twinx()
ln2 = ax2b.plot(t, S_P, color=color_P, linewidth=2, linestyle='--', label='External P (μmol/L)')
ln3 = ax2b.plot(t, Q / Q_min, color=color_Q, linewidth=2, linestyle=':', label='Q/Q_min (cell quota ratio)')
ax2b.set_ylabel('P concentration / Quota ratio', color='white', fontsize=11)
ax2b.tick_params(colors='white')
for sp in ax2b.spines.values(): sp.set_color('#2dd4bf')

lns = ln1 + ln2 + ln3
labs = [l.get_label() for l in lns]
ax2.legend(lns, labs, fontsize=8, facecolor='#1a1a2e', edgecolor='#2dd4bf', labelcolor='white', loc='center right')
ax2.set_xlabel('Time (days)', color='white', fontsize=11)
ax2.set_title('Droop Model: P Limitation with Luxury Uptake\nμ = μ_max(1 - Q_min/Q)', color='white', fontsize=12, fontweight='bold')
ax2.set_facecolor('#0a0a1a')
ax2.tick_params(colors='white')
ax2.grid(True, alpha=0.15, color='#2dd4bf')
for sp in ax2.spines.values(): sp.set_color('#2dd4bf')

# Plot 3: Wetland methane balance
ax3 = axes[1, 0]
ax3.fill_between(T_wetland, 0, P_CH4, alpha=0.2, color='#f87171', label='Gross methanogenesis')
ax3.fill_between(T_wetland, 0, methanotrophy, alpha=0.2, color='#2dd4bf', label='Methanotrophy (oxidation)')
ax3.plot(T_wetland, wt_shallow, color='#fbbf24', linewidth=2, linestyle='--', label='Net: shallow water table')
ax3.plot(T_wetland, net_CH4, color='#f87171', linewidth=2.5, label='Net: normal water table')
ax3.plot(T_wetland, wt_deep, color='#60a5fa', linewidth=2, linestyle='--', label='Net: deep water table')
ax3.axhline(0, color='gray', linewidth=0.5)
ax3.set_xlabel('Temperature (°C)', color='white', fontsize=11)
ax3.set_ylabel('CH₄ flux (mg CH₄/m²/day)', color='white', fontsize=11)
ax3.set_title('Wetland Methane Balance\nProduction (Q₁₀≈4) vs Consumption', color='white', fontsize=12, fontweight='bold')
ax3.legend(fontsize=8, facecolor='#1a1a2e', edgecolor='#2dd4bf', labelcolor='white')
ax3.set_facecolor('#0a0a1a')
ax3.tick_params(colors='white')
ax3.grid(True, alpha=0.15, color='#2dd4bf')
for sp in ax3.spines.values(): sp.set_color('#2dd4bf')

# Plot 4: Gibbs free energy cascade
ax4 = axes[1, 1]
colors_bar = ['#2dd4bf', '#60a5fa', '#fbbf24', '#f87171']
bars = ax4.bar(steps, [-dg for dg in dG_values], color=colors_bar, alpha=0.85, edgecolor='white', linewidth=0.8)
for bar, dg, metal in zip(bars, dG_values, metals):
    ax4.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 5,
             f'{dg} kJ/mol', ha='center', va='bottom', color='white', fontsize=10, fontweight='bold')
    ax4.text(bar.get_x() + bar.get_width()/2, bar.get_height()/2,
             metal, ha='center', va='center', color='#0a0a1a', fontsize=9, fontweight='bold')

ax4.set_ylabel('-ΔG°\' (kJ/mol)', color='white', fontsize=11)
ax4.set_title('Denitrification Energy Cascade\nMetalloenzyme at Each Step', color='white', fontsize=12, fontweight='bold')
ax4.set_facecolor('#0a0a1a')
ax4.tick_params(colors='white')
ax4.grid(True, alpha=0.15, color='#2dd4bf', axis='y')
for sp in ax4.spines.values(): sp.set_color('#2dd4bf')

# Add cumulative line on twin axis
ax4b = ax4.twinx()
positions = np.arange(len(steps))
ax4b.plot(positions, cumulative[1:], 'D-', color='#a78bfa', markersize=8, linewidth=2, label='Cumulative ΔG°\'')
ax4b.set_ylabel('Cumulative ΔG°\' (kJ/mol)', color='#a78bfa', fontsize=11)
ax4b.tick_params(colors='white')
ax4b.legend(fontsize=9, facecolor='#1a1a2e', edgecolor='#a78bfa', labelcolor='white')
for sp in ax4b.spines.values(): sp.set_color('#2dd4bf')

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

# ── Numerical Output ──
print("=== Advanced Soil Biogeochemistry Simulations ===\n")

print("--- Nitrogen Cascade (Denitrification) ---")
peak_idx = np.argmax(N2O_total)
print(f"Peak N₂O emission at WFPS = {WFPS[peak_idx]*100:.0f}%")
print(f"At WFPS < 60%: nitrification dominates N₂O production")
print(f"At WFPS 60-80%: denitrification with incomplete reduction (high N₂O)")
print(f"At WFPS > 80%: complete denitrification to N₂ (low N₂O)\n")

for step, dg in zip(steps, dG_values):
    step_clean = step.replace('\n', ' ')
    print(f"  {step_clean}: ΔG°' = {dg} kJ/mol")
print(f"  Total cascade: ΔG°' = {sum(dG_values)} kJ/mol\n")

print("--- Droop Model (P Limitation) ---")
print(f"Initial P: {S_P[0]:.1f} μmol/L")
print(f"Final P: {S_P[-1]:.3f} μmol/L")
print(f"Peak biomass: {B.max()/1e4:.0f} × 10⁴ cells/mL")
print(f"Final Q/Q_min: {Q[-1]/Q_min:.2f} (growth ceases at Q/Q_min → 1)\n")

print("--- Wetland Methane Balance ---")
print(f"At 15°C: gross production = {P_ref:.0f} mg CH₄/m²/day")
print(f"At 25°C: gross production = {P_ref * Q10_meth ** 1:.0f} mg CH₄/m²/day")
print(f"At 35°C: gross production = {P_ref * Q10_meth ** 2:.0f} mg CH₄/m²/day")
print(f"Methanotrophy oxidizes ~{f_ox_base*100:.0f}% of gross production")
print(f"Q₁₀(methanogenesis) ≈ {Q10_meth:.1f} — highly temperature-sensitive")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

References

Parton, W.J. et al. (1987). Analysis of factors controlling soil organic matter levels in Great Plains grasslands. Soil Science Society of America Journal, 51(5), 1173-1179.
Schimel, J.P. & Schaeffer, S.M. (2012). Microbial control over carbon cycling in soil. Frontiers in Microbiology, 3, 348.
Hoffman, B.M. et al. (2014). Mechanism of nitrogen fixation by nitrogenase: the next stage. Chemical Reviews, 114(8), 4041-4062.
Allison, S.D. et al. (2010). Soil-carbon response to warming dependent on microbial physiology. Nature Geoscience, 3(5), 336-340.
Fierer, N. (2017). Embracing the unknown: disentangling the complexities of the soil microbiome. Nature Reviews Microbiology, 15(10), 579-590.
Louca, S. et al. (2018). Function and functional redundancy in microbial systems. Nature Ecology & Evolution, 2(6), 936-943.
Oldroyd, G.E.D. (2013). Speak, friend, and enter: signalling systems that promote beneficial symbiotic associations in plants. Nature Reviews Microbiology, 11(4), 252-263.
Burns, R.G. et al. (2013). Soil enzymes in a changing environment: current knowledge and future directions. Soil Biology and Biochemistry, 58, 216-234.
Firestone, M.K. & Davidson, E.A. (1989). Microbiological basis of NO and N\(_2\)O production and consumption in soil. Exchange of Trace Gases between Terrestrial Ecosystems and the Atmosphere, 47, 7-21.
Schlesinger, W.H. (1997). Biogeochemistry: An Analysis of Global Change. Academic Press.
Thauer, R.K. et al. (2008). Methanogenic archaea: ecologically relevant differences in energy conservation. Nature Reviews Microbiology, 6(8), 579-591.
Redfield, A.C. (1958). The biological control of chemical factors in the environment. American Scientist, 46(3), 230A-221.

Share:X Reddit LinkedIn

←M2: Animal-Plant Interactions M4: Climate-Dependent Biochemistry→