Quantum Mechanics in Biology

import numpy as np
import matplotlib.pyplot as plt

# Physical constants
hbar = 1.055e-34  # J·s
m_p = 1.673e-27   # proton mass (kg)
eV = 1.602e-19    # electron volt in Joules

# Barrier parameters (typical for enzyme catalysis)
V0 = 0.5 * eV     # barrier height (0.5 eV)
a = 1.0e-10       # barrier width (1 Angstrom)

def tunneling_probability(E, V0, a, m):
    """
    Calculate quantum tunneling probability using WKB approximation

    Parameters:
    E  : particle energy (J)
    V0 : barrier height (J)
    a  : barrier width (m)
    m  : particle mass (kg)
    """
    if E >= V0:
        # Classical case: particle goes over barrier
        return 1.0

    # Quantum tunneling probability
    kappa = np.sqrt(2 * m * (V0 - E)) / hbar
    T = np.exp(-2 * kappa * a)

    return T

# Energy range (from 0 to barrier height)
E_range = np.linspace(0.01 * eV, V0, 100)

# Calculate tunneling probabilities
T_proton = [tunneling_probability(E, V0, a, m_p) for E in E_range]

# Compare with classical expectation (step function)
T_classical = [1.0 if E >= V0 else 0.0 for E in E_range]

# Calculate specific examples
print("Quantum Tunneling in Enzyme Catalysis")
print("=" * 60)
print(f"Barrier height: {V0/eV:.2f} eV")
print(f"Barrier width: {a*1e10:.2f} Å")
print(f"Particle: proton (mass = {m_p:.2e} kg)")
print()

# Test different energies
test_energies = [0.1, 0.2, 0.3, 0.4, 0.45]
print("Tunneling Probabilities:")
print("-" * 60)
for e_frac in test_energies:
    E = e_frac * V0
    T = tunneling_probability(E, V0, a, m_p)
    print(f"E = {e_frac:.2f}V₀ = {E/eV:.3f} eV:")
    print(f"  T_quantum = {T:.2e} ({T*100:.4f}%)")
    print(f"  T_classical = 0 (forbidden)")
    enhancement = T / 1e-100 if T > 0 else 0
    print(f"  Enhancement: {T:.2e} (quantum enables the reaction!)")
    print()

# Temperature effect on reaction rates
print("Temperature Dependence:")
print("-" * 60)
kB = 1.381e-23  # Boltzmann constant
temps = [273, 300, 310, 350]  # Kelvin
for T_kelvin in temps:
    # Thermal energy
    E_thermal = kB * T_kelvin
    # Tunneling probability at thermal energy
    if E_thermal < V0:
        T_tunnel = tunneling_probability(E_thermal, V0, a, m_p)
        print(f"T = {T_kelvin}K: E_thermal = {E_thermal/eV:.4f} eV")
        print(f"  Tunneling prob: {T_tunnel:.2e}")
    print()

# Visualization
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# Plot 1: Tunneling probability vs energy
ax1.semilogy(E_range/eV, T_proton, 'b-', linewidth=2, label='Quantum (proton)')
ax1.plot(E_range/eV, T_classical, 'r--', linewidth=2, label='Classical')
ax1.axvline(x=V0/eV, color='gray', linestyle=':', alpha=0.5, label='Barrier height')
ax1.set_xlabel('Energy (eV)')
ax1.set_ylabel('Transmission Probability')
ax1.set_title('Quantum Tunneling vs Classical Barrier')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Plot 2: Barrier shape with wave function
x = np.linspace(-2*a, 3*a, 1000)
V = np.where((x >= 0) & (x <= a), V0, 0)

ax2.plot(x*1e10, V/eV, 'k-', linewidth=2, label='Potential barrier')
ax2.axhline(y=0.3, color='blue', linestyle='--', alpha=0.7, label='E = 0.3V₀')
ax2.fill_between(x*1e10, 0, V0/eV, where=((x>=0) & (x<=a)),
                  alpha=0.2, color='red', label='Classically forbidden')
ax2.set_xlabel('Position (Å)')
ax2.set_ylabel('Energy (eV)')
ax2.set_title('Energy Barrier in Biological System')
ax2.legend()
ax2.grid(True, alpha=0.3)
ax2.set_ylim(-0.1, 0.7)

plt.tight_layout()
plt.savefig('quantum_tunneling_biology.png', dpi=150, bbox_inches='tight')
print("Plot saved as 'quantum_tunneling_biology.png'")