STRC Research

Recipe — LRA Method for pKa Shift Calculation

5 min read

Recipe — LRA Method for pKa Shift Calculation

P1 recipe — Eq. (12.36) and (12.40)–(12.41) from 2007-chipot-free-energy-calculations-book §12.3.2 (Simonson). The Linear Response Approximation (LRA) method computes the protein-induced pKa shift of a titratable residue from two MD simulations — one of the protonated state, one of the deprotonated state — without running an alchemical λ scan. Use this whenever a charge-state change in a protein needs a defensible free-energy estimate. Directly applicable to the STRC E1659A mutation (charge change), to validating the STRC Electrostatic Analysis E1659A ΔΔG=8.4 kcal/mol number, and to any future titration question in stereocilin.

The core LRA equation (Eq. 12.36, verbatim)

For a single perturbing point charge q inserted onto a side-chain atom:

ΔA = (q / 2) · ( ⟨V⟩_reactant + ⟨V⟩_product )                 [Eq. 12.36]

For a realistic proton modeled as a set of incremental charge shifts {Δq_i} on the side-chain atoms (Eq. 12.40, the general form to use):

ΔA = (1/2) ∑_{i ∈ A} Δq_i · ( ⟨V_i^{reactant}⟩_0 + ⟨V_i^{product}⟩_1 )

where:

  • ⟨V_i⟩_0 = ensemble-average electrostatic potential on side-chain atom i, taken from MD of the reactant (unprotonated) state.
  • ⟨V_i⟩_1 = ensemble-average electrostatic potential on side-chain atom i, taken from MD of the product (protonated) state.
  • Δq_i = the partial-charge increment on atom i between the two protonation states (taken from the force field).
  • A = the side chain undergoing protonation.

— [Chipot & Pohorille 2007, p. 411, §12.3.2 Eq. 12.36 and 12.40]

The pKa shift relative to a model compound (e.g., methylimidazole for histidine, methylglutamate for glutamate) is:

pKa_protein − pKa_model = (1 / (2.303 · k_B · T)) · (ΔA_protein − ΔA_model)        [Eq. 12.41]

— [Chipot & Pohorille 2007, p. 411, Eq. 12.41]

Why LRA is the right tool

  • Quadratic free-energy curve assumption: the book derives this from the dielectric response being linear in the perturbing charge (Eq. 12.21–12.32). For modest charge perturbations (single-residue ionization) this is well-tested.
  • Two MDs replace a full λ scan: instead of 10–20 alchemical windows you need exactly two equilibrium MDs (one per protonation state) plus electrostatic-potential post-processing.
  • Midpoint shortcut (Eq. 12.38): if you can run a third MD at fractional charge (λ = 0.5), then ΔA = ΔA_static^{midpoint} — one simulation gives ΔA. This is the single-cheapest pKa estimator in the book.
  • Validates against MDFE: Simonson reports LRA gives near-MDFE quality on shifted pKa (within ~1 pKa unit on most proteins, error of 3 pKa units on one outlier in thioredoxin §12.6.4). The book pairs it with PB/LRA where the explicit-solvent step is replaced by Poisson–Boltzmann; comparable accuracy at lower cost.

Protocol

  1. Pick the model compound — for E1659A in STRC, the natural model is acetic acid or N-methylacetamide for the glutamate side chain. Run a 50–100 ns MD of the model in TIP3P water; record the average electrostatic potential at the carboxylate oxygens.
  2. Run reactant MD of the protein in the unprotonated state (charged Asp/Glu, or unprotonated His). 100+ ns; converged side-chain conformation; standard force field (AMBER FF19SB / CHARMM36m).
  3. Run product MD of the protein in the protonated state (neutral Asp/Glu COOH, or protonated His⁺). Same length and force field.
  4. Compute ⟨V_i⟩_0 and ⟨V_i⟩_1 at every titrating-side-chain atom, averaged over each trajectory. Tools: VMD volmap, GROMACS gmx potential, or a custom NumPy script over prmtop + nc frames.
  5. Plug into Eq. 12.40 with the force-field charge increments Δq_i.
  6. Repeat steps 2–5 for the model compound in solution.
  7. Compute ΔΔA = ΔA_protein − ΔA_model and convert to pKa shift via Eq. 12.41.

Quality checks (book §12.6.4)

  • Charge set vs dielectric consistency — if you also run a continuum (PB) calculation as a cross-check, use solute dielectric ε ∈ [1, 2]. ε = 4 amounts to double-counting protein relaxation already captured in the MD (2007-chipot-free-energy-calculations-book §12.6.4, p. 461).
  • Variance(η)/2 cross-check (Eq. 12.32, fluctuation–dissipation theorem) — the relaxation free energy can be computed two ways: as −β/2 · Var(η) and as (⟨η⟩_1 − ⟨η⟩_0)/2. They must agree if the linear-response assumption holds. If they disagree by >5 kcal/mol the system has nonlinear protonation-induced relaxation and LRA is suspect — switch to alchemical λ.
  • Endpoint structure check — book §12.6.4 emphasises that for highly relaxing systems both endpoint structures must be sampled; using only the X-ray structure can give 4–5 kcal/mol errors per [Archontis et al.], cited at p. 459.

How to use this in STRC

  • STRC Electrostatic Analysis E1659A: the existing 8.4 kcal/mol multi-method ΔΔG number can be cross-validated via LRA at a fraction of the cost of a full alchemical scan. Proposal: run 100 ns MD of WT and of E1659A in TIP3P water with AMBER FF19SB (already the h01 force field per pharmacochaperone); apply Eq. 12.40; compare to existing 8.4 kcal/mol.
  • h01 pharmacochaperone screening — when a candidate ligand changes the protonation environment of a nearby residue (e.g., binds near K1141 or D1657), LRA gives the per-ligand pKa shift from the same MD trajectories used for MM-PBSA scoring — no extra simulations.
  • h26 disulfide cysteines — each Cys gain shifts the local pKa of nearby His/Glu. LRA gives the shift per design directly from the AF3-relaxed MDs.

Limitations

  • Strictly linear-response: large protonation-induced reorganisation (book’s thioredoxin Asp26 case, ΔA_rlx = −56 kcal/mol) is the failure mode. Always run the variance cross-check.
  • Force-field charge set must match: Δq_i must come from the same FF used for the MD. Switching from AMBER to CHARMM mid-pipeline invalidates the calculation.
  • Cannot capture proton-transfer kinetics — purely a thermodynamic shift.

Connections

Graph View

Backlinks