Force Field¶

The AdsPro force field has four components:

Morse potential — protein residue ↔ NP surface bead interaction
Debye-Hückel potential — long-range electrostatics (residue ↔ surface bead)
Gō model — protein internal energy (native contacts)
Harmonic backbone restraints — chain connectivity

Morse Potential: Protein–Surface Interaction¶

The short-range non-electrostatic interaction between protein residue i and surface bead j is modeled by a Morse potential, which provides a physically motivated form: steep repulsion at short range, a well-defined attractive minimum, and smooth asymptotic decay to zero at large distance.

Attractive Branch (ε < 0: residue binds surface)¶

\[ V_\text{Morse}(r) = \varepsilon \left[1 - e^{-\alpha(r - r_\text{min})}\right]^2 - \varepsilon \]

At the minimum distance \(r = r_\text{min}\): \(V = -\varepsilon\) (well depth).
As \(r \to \infty\): \(V \to 0\).
The "repulsive wall" at \(r < r_\text{min}\) rises steeply as \(e^{-2\alpha(r-r_\text{min})}\).

Repulsive Branch (ε > 0: residue is repelled)¶

For residues repelled by the surface (e.g., ASP/GLU on negatively charged silica), a purely repulsive exponential is used:

\[ V_\text{repulsive}(r) = \varepsilon \cdot e^{-2\alpha(r - r_\text{min})} \]

This decays to zero as the residue moves away from the surface, without creating a spurious attractive well.

Parameters¶

Parameter	Symbol	Value	Source
Stiffness	α	8.0 nm⁻¹	Fit to amino acid adsorption data
Minimum distance	r_min	0.0 nm	Hard contact
Cutoff	r_cut	1.5 nm	Beyond this V = 0
Well depth	ε	residue-specific	Bag et al. (2020) chromatography

Parameterization of ε from Chromatography¶

The Morse well depth ε per residue type is derived from the chromatography experiments of Bag et al. (2020) (ChemPhysChem 21:2347). Individual amino acids were loaded onto a Q3 silica HPLC column at nanomolar to low micromolar concentration (Henry's law regime) and eluted with physiological buffer (50 mM MOPS, pH 7.4, 150 mM NaCl).

The measured quantity is the Henry coefficient K_H — the equilibrium partition coefficient between surface and solution at trace concentration:

\[ K_H = \frac{V_\text{net}}{V_\text{dead}} \]

This directly yields the residue-level free energy of adsorption:

\[ \varepsilon = \Delta G^\circ_\text{ads} = -RT \ln K_H \]

Values for the 20 amino acids at the Q3 silica surface:

Residue	ε (kJ/mol)	Mechanism
ARG	−15	Bidentate guanidinium–SiO⁻ ion pairing
LYS	−9	NH₃⁺–SiO⁻ electrostatics + H-bond
HIS	−1.5	Partial protonation (~12% at pH 7.4)
ASP	+3	COO⁻ repelled by SiO⁻
GLU	+3	COO⁻ repelled by SiO⁻
All others	0	No detectable interaction

Physical insight

Silica adsorption is driven by charged residues only. Hydrophobic residues (TRP, PHE, ILE, LEU, VAL) show zero interaction with silica at physiological pH because the silica surface is polar and hydrophilic. This is fundamentally different from hydrophobic NP surfaces — see Extending to New Surfaces.

Henry's law regime

The experiment must be conducted at nanomolar to low micromolar amino acid concentrations to ensure K_H is concentration-independent. At higher concentrations, surface sites become saturated, K_H decreases, and the measured ΔG is no longer the single-molecule adsorption energy needed by AdsPro.

Debye-Hückel Electrostatics: Residue–Surface Interaction¶

Long-range electrostatic interactions between charged protein residues and charged surface beads are computed via the screened Coulomb (Debye-Hückel) potential:

\[ V_\text{DH}(r) = \frac{q_i \, q_j \, K_\text{Coulomb}}{\varepsilon_r \, r} \, e^{-r/\lambda_D} \]

where: - \(K_\text{Coulomb} = 138.935\) kJ·nm/(mol·e²) — Coulomb constant in MD units - \(\varepsilon_r = 78.4\) — relative permittivity of water - \(\lambda_D\) — Debye screening length (computed from ionic strength and temperature) - Cutoff: \(3\lambda_D\) — beyond this, the interaction is negligibly small

The exponential factor \(e^{-r/\lambda_D}\) encodes the screening by the ionic atmosphere. At distance \(r = \lambda_D\), the interaction is reduced to \(1/e \approx 37\%\) of its unscreened value.

Residue–Residue Electrostatics¶

The same Debye-Hückel potential governs electrostatics between charged protein residues (e.g., ARG 14 repelling LYS 21 within the protein, or ARG 14 attracting GLU 35). These intra-protein electrostatic interactions are part of the internal energy that changes upon adsorption.

Gō Model: Protein Internal Energy¶

Native contacts (detected from the PDB, see Protein Model) are assigned a 12-6 Lennard-Jones attractive potential:

\[ V_{G\bar{o}}(r) = \varepsilon_{G\bar{o}} \left[ \left(\frac{r_0}{r}\right)^{12} - 2\left(\frac{r_0}{r}\right)^{6} \right] \]

Minimum at \(r = r_0\) (native distance): \(V = -\varepsilon_{G\bar{o}}\)
Repulsive wall at \(r < r_0\): \(V \to +\infty\)
Attractive tail at \(r > r_0\): \(V \to 0\)

Non-native pairs interact through a WCA (Weeks-Chandler-Andersen) purely repulsive potential, which is the repulsive branch of a 12-6 LJ potential:

\[ V_\text{WCA}(r) = \begin{cases} 4\varepsilon \left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\right] + \varepsilon & r < 2^{1/6}\sigma \\ 0 & r \geq 2^{1/6}\sigma \end{cases} \]

This ensures excluded volume (no two beads can occupy the same space) without any spurious attraction between non-native pairs.

Harmonic Backbone Restraints¶

Adjacent Cα–Cα+1 pairs are held at their native PDB bond lengths by harmonic spring restraints:

\[ V_\text{backbone}(r) = \frac{1}{2} K_\text{bond} (r - r_0)^2 \]

with \(K_\text{bond} = 8000\) kJ/(mol·nm²).

This stiffness corresponds to a force of 160 kJ/mol·nm per 0.02 nm (2%) elongation — strong enough to prevent unphysical chain stretching while allowing thermal fluctuations of ≈ 0.02 nm (√(kBT/K) at 310 K).

Not included in energy balance

Backbone spring energy is deliberately excluded from the reported energy summary. These springs are constraint-like forces that maintain chain topology — they are not physical interactions contributing to ΔG_ads. Only Morse, Debye-Hückel, and Gō energies appear in the energy balance.

Force Cutoffs and Performance¶

Interaction	Cutoff	Justification
Morse	1.5 nm	Exactly zero beyond this by construction
Debye-Hückel	3 λ_D	At 3λ_D, interaction is exp(−3) ≈ 5% of contact value
Gō contacts	Native distance × 3	LJ is negligible at 3r_0

The active zone optimization computes which surface beads are within cutoff of any protein residue at each step, avoiding O(N_residues × N_surface_beads) loops. For a 71.3 nm NP (~64,000 beads) with a 129-residue protein, this reduces the per-step surface interaction computation from ~8 million pairs to a few hundred.