Computational Biophysics

Derivation: The Verlet Algorithm

Start with Taylor expansions of the position forward and backward in time:

$$\mathbf{r}(t + \Delta t) = \mathbf{r}(t) + \mathbf{v}(t)\Delta t + \frac{1}{2}\mathbf{a}(t)\Delta t^2 + \frac{1}{6}\mathbf{b}(t)\Delta t^3 + O(\Delta t^4)$$

$$\mathbf{r}(t - \Delta t) = \mathbf{r}(t) - \mathbf{v}(t)\Delta t + \frac{1}{2}\mathbf{a}(t)\Delta t^2 - \frac{1}{6}\mathbf{b}(t)\Delta t^3 + O(\Delta t^4)$$

Adding these two equations, the odd-order terms cancel:

$$\boxed{\mathbf{r}(t + \Delta t) = 2\mathbf{r}(t) - \mathbf{r}(t - \Delta t) + \mathbf{a}(t)\Delta t^2 + O(\Delta t^4)}$$

where $\mathbf{a}(t) = \mathbf{F}(t)/m$. This is the Verlet algorithm. Key properties:

• Time-reversible (symmetric under $t \to -t$)
• Symplectic (preserves phase-space volume, no energy drift)
• Local error $O(\Delta t^4)$ for positions
• Does not explicitly use velocities (they can be computed from finite differences)

Velocities, needed for kinetic energy and temperature, are obtained as:

$$\mathbf{v}(t) = \frac{\mathbf{r}(t + \Delta t) - \mathbf{r}(t - \Delta t)}{2\Delta t} + O(\Delta t^2)$$

Derivation: Velocity Verlet Algorithm

A more practical variant that maintains positions and velocities at the same time step:

$$\mathbf{r}(t + \Delta t) = \mathbf{r}(t) + \mathbf{v}(t)\Delta t + \frac{1}{2}\mathbf{a}(t)\Delta t^2$$

$$\mathbf{v}(t + \Delta t) = \mathbf{v}(t) + \frac{1}{2}[\mathbf{a}(t) + \mathbf{a}(t + \Delta t)]\Delta t$$

This is mathematically equivalent to the Verlet algorithm but avoids numerical precision issues. The implementation proceeds in three steps: (1) compute $\mathbf{r}(t+\Delta t)$ and half-step velocity $\mathbf{v}(t+\Delta t/2) = \mathbf{v}(t) + \frac{1}{2}\mathbf{a}(t)\Delta t$; (2) compute forces at the new positions to get $\mathbf{a}(t+\Delta t)$; (3) complete velocity update $\mathbf{v}(t+\Delta t) = \mathbf{v}(t+\Delta t/2) + \frac{1}{2}\mathbf{a}(t+\Delta t)\Delta t$.

Derivation: Time Step Constraint

The integration is stable only if the time step $\Delta t$ is much smaller than the shortest oscillation period in the system. For a harmonic oscillator with frequency$\omega$, stability requires:

$$\Delta t < \frac{2}{\omega_{\max}}$$

In biomolecular systems, the highest-frequency motions are C–H, N–H, and O–H bond stretches with $\omega \approx 3000$ cm$^{-1}$($\sim 10^{14}$ s$^{-1}$), giving a period of$\sim 10$ fs. A safe time step is:

$$\boxed{\Delta t \lesssim \frac{1}{\omega_{\max}} \approx 1 \text{ fs}}$$

SHAKE/LINCS constraint algorithms: By constraining bonds to hydrogen atoms (fixing their length), these fast vibrations are removed, allowing a 2 fs time step — doubling the effective sampling rate. Hydrogen mass repartitioning (HMR) further extends this to 4 fs by redistributing mass from heavy atoms to hydrogens, slowing the constrained vibrations.

Computational cost: A typical explicit-solvent MD system has $\sim 10^5$ atoms. Each time step requires $O(N \log N)$ operations (using PME for electrostatics). At 2 fs/step, simulating 1 $\mu$s requires$5 \times 10^8$ steps. On modern GPUs, this takes about 1 day for a 100k-atom system.

Derivation: The Molecular Mechanics Potential Energy Function

The total potential energy is decomposed into bonded (intramolecular) and nonbonded terms:

$$V_{\text{total}} = V_{\text{bonds}} + V_{\text{angles}} + V_{\text{dihedrals}} + V_{\text{vdW}} + V_{\text{elec}}$$

1. Bond stretching (harmonic):

$$V_{\text{bonds}} = \sum_{\text{bonds}} K_b(b - b_0)^2$$

where $b$ is the bond length, $b_0$ the equilibrium value, and$K_b$ the force constant (typically $200\text{--}500$ kcal/mol/Å$^2$). The force is $F = -dV/db = -2K_b(b-b_0)$, a restoring force toward equilibrium.

2. Angle bending (harmonic):

$$V_{\text{angles}} = \sum_{\text{angles}} K_\theta(\theta - \theta_0)^2$$

with force constants typically $30\text{--}100$ kcal/mol/rad$^2$.

3. Dihedral (torsional) potentials:

$$V_{\text{dihedrals}} = \sum_{\text{dihedrals}} K_\phi[1 + \cos(n\phi - \delta)]$$

where $\phi$ is the dihedral angle, $n$ is the periodicity (typically 1–6), $\delta$ is the phase, and $K_\phi$ the barrier height. This term governs rotational barriers (e.g., Ramachandran $\phi/\psi$ angles in proteins).

4. van der Waals interactions (Lennard-Jones):

$$V_{\text{vdW}} = \sum_{i<j} 4\varepsilon_{ij}\left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12} - \left(\frac{\sigma_{ij}}{r_{ij}}\right)^6\right]$$

The $r^{-12}$ term models Pauli repulsion; the $r^{-6}$ term captures London dispersion forces. Combination rules: $\sigma_{ij} = (\sigma_i + \sigma_j)/2$,$\varepsilon_{ij} = \sqrt{\varepsilon_i \varepsilon_j}$ (Lorentz-Berthelot).

5. Electrostatic interactions (Coulomb):

$$V_{\text{elec}} = \sum_{i<j} \frac{q_i q_j}{4\pi\varepsilon_0 r_{ij}}$$

Partial atomic charges $q_i$ are derived from quantum mechanical calculations (e.g., RESP fitting: the electrostatic potential from an ab initio calculation is reproduced by atom-centered point charges, constrained to reproduce dipole moments and minimize non-physical charge magnitudes).

Particle Mesh Ewald (PME): Long-range electrostatics in periodic systems are handled by splitting the Coulomb sum into short-range (real space,$O(N)$) and long-range (reciprocal space, FFT, $O(N\log N)$) contributions, reducing the cost from $O(N^2)$ to $O(N\log N)$.

Derivation: Thermodynamic Integration (TI)

Consider transforming a system from state A (Hamiltonian $H_A$) to state B (Hamiltonian $H_B$) via a coupling parameter $\lambda$:

$$H(\lambda) = (1-\lambda)H_A + \lambda H_B$$

The free energy of state $\lambda$ is $G(\lambda) = -k_BT\ln Z(\lambda)$. Taking the derivative with respect to $\lambda$:

$$\frac{dG}{d\lambda} = \frac{1}{Z(\lambda)}\int \frac{\partial H}{\partial\lambda} e^{-H(\lambda)/k_BT} d\mathbf{r}\, d\mathbf{p} = \left\langle\frac{\partial H}{\partial\lambda}\right\rangle_\lambda$$

Integrating from $\lambda = 0$ to $\lambda = 1$:

$$\boxed{\Delta G = G_B - G_A = \int_0^1 \left\langle\frac{\partial H}{\partial\lambda}\right\rangle_\lambda d\lambda}$$

In practice, simulations are run at discrete $\lambda$ values (typically 11–21 windows), the ensemble average $\langle\partial H/\partial\lambda\rangle_\lambda$ is computed at each, and the integral is evaluated numerically (trapezoidal or Gaussian quadrature).

Derivation: Free Energy Perturbation (FEP)

The Zwanzig formula (1954) relates the free energy difference to a single-ensemble average:

$$\Delta G = G_B - G_A = -k_BT \ln\left\langle\exp\left(-\frac{H_B - H_A}{k_BT}\right)\right\rangle_A$$

Proof: Starting from the ratio of partition functions:

$$e^{-\Delta G/k_BT} = \frac{Z_B}{Z_A} = \frac{\int e^{-H_B/k_BT}d\Gamma}{\int e^{-H_A/k_BT}d\Gamma} = \frac{\int e^{-(H_B-H_A)/k_BT} e^{-H_A/k_BT}d\Gamma}{\int e^{-H_A/k_BT}d\Gamma}$$

$$= \left\langle e^{-(H_B-H_A)/k_BT}\right\rangle_A$$

Taking the logarithm gives the Zwanzig formula. The exponential averaging makes FEP sensitive to rare configurations where $H_B - H_A$ is large and negative (poor overlap between ensembles). This limits its use to small perturbations ($|\Delta G| \lesssim 2\text{--}3 k_BT$) unless staged.

Umbrella Sampling and WHAM

For computing the potential of mean force (PMF) along a reaction coordinate $\xi$(e.g., ligand unbinding distance), umbrella sampling applies harmonic biasing potentials to ensure adequate sampling at each point:

$$w_i(\xi) = \frac{1}{2}k(\xi - \xi_i)^2$$

The biased probability distributions are combined using the Weighted Histogram Analysis Method (WHAM) to reconstruct the unbiased PMF:

$$P_{\text{unbiased}}(\xi) = \frac{\sum_i n_i P_i^{\text{biased}}(\xi) e^{w_i(\xi)/k_BT}}{\sum_i n_i e^{-(f_i - w_i(\xi))/k_BT}}$$

where $f_i$ are free energies of each window, determined self-consistently. The PMF is then $W(\xi) = -k_BT \ln P_{\text{unbiased}}(\xi) + \text{const}$.

Thermodynamic Cycles for Binding Free Energies

Computing absolute binding free energies directly is challenging because of the large conformational entropy change upon binding. Instead, relative binding free energies are computed using a thermodynamic cycle. For two ligands A and B binding to the same protein P:

$$\Delta\Delta G_{\text{bind}} = \Delta G_{\text{bind}}^B - \Delta G_{\text{bind}}^A = \Delta G_{\text{complex}}^{A\to B} - \Delta G_{\text{solvent}}^{A\to B}$$

The alchemical transformations $A \to B$ in complex and in solvent are computed via TI or FEP with $\lambda$ windows. The double free energy difference$\Delta\Delta G$ benefits from cancellation of systematic errors (force field inaccuracies, finite-size effects) that are similar in both legs of the cycle.

Bennett Acceptance Ratio (BAR): An optimal estimator that uses samples from both forward ($A \to B$) and backward ($B \to A$) perturbations:

$$\Delta G = k_BT\ln\frac{\langle f(\Delta H + C)\rangle_B}{\langle f(-\Delta H - C)\rangle_A} + C$$

where $f(x) = 1/(1 + e^{x/k_BT})$ is the Fermi function and $C$ is determined self-consistently. BAR achieves the lowest variance among all estimators using data from two states.

Practical Considerations in Free Energy Calculations

Soft-core potentials: When transforming atoms (appearing/disappearing), the $r^{-12}$ repulsion causes numerical instabilities at small $\lambda$. Soft-core Lennard-Jones potentials replace $r$ with$(\alpha\sigma^6\lambda^p + r^6)^{1/6}$ to avoid singularities.

Long-range corrections: Finite-size effects from periodic boundary conditions require corrections for charged species. The Rocklin correction accounts for the interaction of the net charge with its periodic images:$\Delta G_{\text{corr}} \propto q^2 / (L\varepsilon)$.

Convergence assessment: Free energy estimates should be monitored as a function of simulation time, and forward/backward estimates should agree within statistical uncertainty. Overlap between adjacent $\lambda$ windows (measured by the overlap of energy distributions) should be $> 10\%$ for reliable estimates.

Derivation: Replica Exchange Molecular Dynamics (REMD)

In REMD (parallel tempering), $M$ replicas of the system are simulated simultaneously at different temperatures $T_1 < T_2 < \cdots < T_M$. Periodically, exchanges of configurations between adjacent temperatures are attempted.

For replicas $i$ and $j$ at temperatures $T_i$ and$T_j$ with potential energies $E_i$ and $E_j$, the joint probability in the canonical ensemble is:

$$P(\mathbf{r}_i, T_i; \mathbf{r}_j, T_j) \propto \exp\left(-\frac{E_i}{k_BT_i} - \frac{E_j}{k_BT_j}\right)$$

After exchange, the probability becomes:

$$P(\mathbf{r}_j, T_i; \mathbf{r}_i, T_j) \propto \exp\left(-\frac{E_j}{k_BT_i} - \frac{E_i}{k_BT_j}\right)$$

The Metropolis acceptance criterion for detailed balance:

$$\boxed{P_{\text{accept}} = \min\left(1, \exp\left(\Delta\beta \cdot \Delta E\right)\right)}$$

where $\Delta\beta = 1/(k_BT_i) - 1/(k_BT_j)$ and$\Delta E = E_i - E_j$. The exchange is favorable when a high-energy configuration moves to a higher temperature (where high energies are more probable) and a low-energy configuration moves to a lower temperature.

Temperature spacing: For adequate exchange rates ($\sim 20\text{--}30\%$ acceptance), temperatures should be spaced geometrically:$T_{i+1}/T_i = \text{const}$. The number of replicas scales as$M \propto \sqrt{N_{\text{atoms}}}$, limiting REMD to systems of$\lesssim 50{,}000$ atoms.

Derivation: Metadynamics

Metadynamics (Laio and Parrinello, 2002) accelerates sampling by progressively adding repulsive Gaussian potentials along chosen collective variables (CVs)$\mathbf{s} = \mathbf{s}(\mathbf{r})$:

$$\boxed{V_{\text{bias}}(\mathbf{s}, t) = \sum_{t' < t} w \exp\left(-\frac{(\mathbf{s} - \mathbf{s}(t'))^2}{2\sigma^2}\right)}$$

where $w$ is the Gaussian height, $\sigma$ the width, and the sum runs over all previously deposited Gaussians. As the simulation progresses, the bias potential fills up free energy minima, eventually enabling the system to escape and explore new regions.

Convergence: In the long-time limit, the bias potential converges to the negative of the free energy surface:

$$V_{\text{bias}}(\mathbf{s}, t \to \infty) \to -F(\mathbf{s}) + C$$

where $F(\mathbf{s})$ is the free energy (PMF) along the CVs and $C$is a constant. This is because once the bias exactly compensates the free energy, the system diffuses freely in CV space, and no more bias is deposited in equilibrium.

Well-tempered metadynamics improves convergence by reducing the Gaussian height over time: $w(t) = w_0 \exp(-V_{\text{bias}}(\mathbf{s}(t), t) / (k_B \Delta T))$, where $\Delta T$ is the bias temperature. The bias converges to$-(\Delta T / (T + \Delta T)) F(\mathbf{s})$.

Derivation: MARTINI Mapping and Speedup

The MARTINI force field (Marrink et al., 2007) maps approximately 4 heavy atoms to 1 CG bead. For example, a phospholipid with $\sim 130$ atoms is represented by$\sim 12$ beads, a reduction factor of $\sim 10$.

Speedup analysis: The computational speedup comes from three factors:

1. Fewer particles: If $N_{\text{CG}} = N_{\text{AA}}/\alpha$(where $\alpha \approx 10$ is the mapping ratio), and the cost of force calculation scales as $O(N)$ per step (with cutoffs), the speedup per step is$\sim \alpha$.

2. Smoother energy landscape: CG potentials are softer, removing fast vibrations. The highest-frequency modes ($\omega_{\max}^{\text{CG}} \approx \omega_{\max}^{\text{AA}} / \sqrt{\alpha}$) allow larger time steps:

$$\Delta t_{\text{CG}} \approx \sqrt{\alpha} \cdot \Delta t_{\text{AA}}$$

For MARTINI, $\Delta t_{\text{CG}} = 20\text{--}30$ fs (vs 2 fs for atomistic).

3. Faster dynamics: CG potentials are smoother, reducing friction and accelerating diffusion by a factor of $\sim 4$ for MARTINI.

The total effective speedup is:

$$\boxed{S_{\text{total}} \approx \alpha \times \sqrt{\alpha} \times f_{\text{dynamics}} \approx 10 \times 3 \times 4 \approx 120}$$

With additional speedup from implicit solvent (removing water beads), the effective speedup can reach $\sim 1000\times$ or more, enabling microsecond-to-millisecond simulations of membrane systems and large assemblies.

Iterative Boltzmann Inversion for CG Potentials

To derive CG interaction potentials that reproduce the correct structural distributions from atomistic simulations, iterative Boltzmann inversion (IBI) starts with an initial guess from the potential of mean force:

$$V_{\text{CG}}^{(0)}(r) = -k_BT \ln g_{\text{target}}(r)$$

where $g_{\text{target}}(r)$ is the radial distribution function from an atomistic simulation of the CG sites. Then iteratively refine:

$$V_{\text{CG}}^{(n+1)}(r) = V_{\text{CG}}^{(n)}(r) + k_BT \ln\frac{g^{(n)}(r)}{g_{\text{target}}(r)}$$

where $g^{(n)}(r)$ is the RDF from a CG simulation using potential$V_{\text{CG}}^{(n)}$. Convergence is typically achieved in 5–20 iterations when $g^{(n)}(r) \to g_{\text{target}}(r)$.

MARTINI Bead Types and Interaction Matrix

The MARTINI force field classifies CG beads into four main types based on chemical character:

• Q (charged): Ions, charged amino acids (Lys, Arg, Glu, Asp)
• P (polar): Water, serine, threonine, backbone amides
• N (nonpolar/intermediate): Mixed character groups
• C (apolar): Hydrophobic groups (alkyl chains, aromatic rings)

Each type has subtypes (e.g., Q0, Qa, Qd, Qda for charge subtypes), giving approximately 18 bead types. The $18 \times 18$ interaction matrix defines LJ parameters ($\varepsilon_{ij}$) between all bead pairs. The interaction levels range from$\varepsilon = 2.0$ kJ/mol (super-repulsive, e.g., charged-apolar) to$\varepsilon = 5.6$ kJ/mol (strong attraction, e.g., polar-polar).

Protein mapping: Each amino acid backbone is represented by one bead (BB), while side chains use 0–4 beads. Gly has no side chain bead; Ala has none (information folded into BB); Trp requires 4 side chain beads to represent the indole ring.

Water in MARTINI: One CG water bead represents 4 atomistic water molecules. To reproduce the dielectric properties of water, the MARTINI model uses a relative dielectric constant of $\varepsilon_r = 15$ (vs 80 for real water), with electrostatic interactions partially screened. Polarizable water models improve this by adding Drude-like particles to each water bead.

Multiscale Simulation Approaches

Backmapping: CG simulations can be converted back to all-atom resolution by inserting atomistic detail into the CG coordinates, followed by energy minimization and short equilibration. This enables a multiscale workflow: (1) CG simulation to explore large-scale conformational space; (2) select representative structures; (3) backmap to all-atom for high-resolution analysis.

Adaptive resolution: Methods like AdResS (Adaptive Resolution Simulation) allow different regions of the simulation box to be treated at different resolutions simultaneously. The active site of a protein can be atomistic while the surrounding solvent is coarse-grained, with a hybrid transition region connecting them.

Machine-learned force fields: Neural network potentials (ANI, NequIP, MACE) trained on quantum mechanical data achieve near-DFT accuracy at a fraction of the cost, enabling simulations with quantum-level accuracy for millions of time steps. These ML potentials are particularly valuable for reactive processes, metalloproteins, and systems where classical force fields are insufficiently accurate.