Plot Interaction Energies - k-ngo/CATMD GitHub Wiki

Plot Interaction Energies

Overview and Methodology

What It Does

This script estimates non-bonded interaction energies between two atom groups over the course of a molecular dynamics trajectory. It computes and visualizes electrostatic, Lennard-Jones (LJ), and hydrophobic contributions using simplified energy models. The tool also attributes average energy contributions to individual residues, enabling insights into which residues drive binding or repel interaction.

This method is intended for qualitative or rapid screening purposes and does not replace rigorous free energy calculations.

How It Works

  • Objective: Identify key energetic contributors and trends in molecular recognition or complex formation.

  • Energy Types:

    • Electrostatic:
      • Based on Coulomb's law with exponential Debye–Hückel screening.
      • Adjustable salt concentration, dielectric constant, and cutoff.
    • Lennard-Jones (LJ):
      • Short-range van der Waals interactions modeled with 12-6 LJ potential.
      • Includes distance cutoff and energy shift correction at cutoff.
    • Hydrophobic:
      • Empirical potential based on residue COM distances and Gaussian-like profile.
      • Applied only between hydrophobic residues.

  • Process:

    • Frame-by-frame calculations for each energy type.
    • Per-residue breakdown for both groups.
    • Rolling average (optional) smooths energy trends over time.
    • Residue contributions filtered based on percent threshold relative to overall distribution.

Energy Models and Equations

Electrostatic Interaction (Debye–Hückel Screening)

$[ E_{\text{elec}} = \frac{332 \cdot q_1 q_2}{\varepsilon \cdot r} \cdot e^{-\kappa r} ]$

  • $( q_1, q_2 )$: partial charges (e)
  • $( r )$: interatomic distance (Å)
  • $( \varepsilon )$: dielectric constant
  • $( \kappa )$: Debye screening constant:

$[ \kappa = 3.289 \cdot \sqrt{c} ]$

where $( c )$ is the salt concentration in mol/L.

  • Default Parameters:
    • k_elec_elec = 332 kcal·Å/mol·e²
    • epsilon_water (default ~5)
    • salt_conc = 0.15 M
    • cutoff_elec = 12 Å

Lennard-Jones Interaction (Shifted 12-6 Potential)

$[ E_{\text{LJ}} = 4\varepsilon \left[\left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^6 \right] - E_{\text{shift}} ]$

  • $( \varepsilon )$: depth of the potential well

  • $( \sigma )$: distance at which potential = 0

  • $( r )$: interatomic distance (Å)

  • $( E_{\text{shift}} )$: energy at cutoff, subtracted to ensure smooth truncation

  • Default Parameters:

    • epsilon_lj = 0.02 kcal/mol
    • sigma_lj = 3.0 Å
    • cutoff_LJ = 12 Å

Hydrophobic Interaction (Gaussian-like Potential)

$[ E_{\text{hp}} = -A \cdot \exp\left( -\frac{(r - r_0)^2}{\delta^2} \right) ]$

  • $( A )$: amplitude (typically tuned to produce ~–30 kcal/mol when summed)

  • $( r_0 )$: optimal interaction distance (Å)

  • $( \delta ):$ interaction width (Å)

  • Applied only to pairs of hydrophobic residues.

  • Default Parameters:

    • A_hp = 23
    • r0_hp = 5.0 Å
    • delta_hp = 1.0 Å
    • cutoff_hp = 12 Å

Configuration and Inputs

Prerequisites

  • Requires a loaded trajectory.

Key Configuration Options

  • Selections:

    • group1_sel, group2_sel: Define receptor and ligand (e.g., protein and toxin).
    • group1_name, group2_name: Used in labeling output plots.

  • Processing:

    • num_threads: Enables parallelization for faster frame analysis.

  • Electrostatics Parameters:

    • salt_conc: Salt concentration (M) for screening factor.
    • epsilon_water: Dielectric constant of the environment.
    • k_elec_elec: Electrostatic constant.
    • cutoff_elec: Distance cutoff for interactions (Å).

  • Lennard-Jones Parameters:

    • epsilon_lj, sigma_lj: Depth and radius for the LJ potential.
    • cutoff_LJ: LJ distance cutoff (Å).

  • Hydrophobic Interaction Parameters:

    • A_hp, r0_hp, delta_hp, cutoff_hp: Amplitude, optimal distance, width, and cutoff for hydrophobic energy.

  • Residue Contribution Filtering:

    • min_contribution_percent: Filters residues below this % of maximum observed contribution.

Output

  • Binding Energy Time-Series Plot:

    • {group1_name}_{group2_name}_binding_energies.png
    • Overlays electrostatic, LJ, hydrophobic, and total (free) energy over time.

  • Per-Residue Energy Contribution Plots:

    • {group1_name}_energy_contributions.png and {group2_name}_energy_contributions.png
    • Horizontal stacked bars grouped by interaction type (electrostatic, LJ, hydrophobic).
    • Only includes residues above specified contribution threshold.

  • Console Output:

    • Average values of each energy term.
    • Debug summaries of per-frame computation (optional).
    • Confirmation of saved plots.

Interpreting the Results

  • Electrostatic:

    • Negative values suggest attractive charge complementarity.
    • Positive spikes may indicate charge repulsion or unscreened buried charges.

  • Lennard-Jones:

    • Usually stabilizing (negative) in the short range.
    • High fluctuations may imply steric clash or contact instability.

  • Hydrophobic:

    • Dominant in buried or dry interfaces.
    • Long-range and smoother compared to LJ.

  • Free Energy:

    • Sum of electrostatic, LJ, and hydrophobic terms.
    • Negative values indicate net favorable binding conditions.

  • Residue Contributions:

    • Identify key stabilizing or destabilizing residues.
    • Guide for mutagenesis, pharmacophore mapping, or interface optimization.

Example Scenarios

Toxin–Channel Binding

  • Scenario: Evaluate toxin binding affinity over simulation.
  • Observation: Hydrophobic contribution increases steadily, while electrostatics remain stable.
  • Interpretation: Suggests buried hydrophobic binding site with stable ionic pairing.

Ligand Entry Pathway

  • Scenario: Ligand gradually enters receptor binding site.
  • Observation: LJ and hydrophobic energy drop over time, electrostatics fluctuate near entry.
  • Interpretation: Consistent with ligand descent into a partially buried pocket.

Hotspot Mapping

  • Scenario: Identify binding hotspots on receptor surface.
  • Observation: Per-residue bar plot shows a few residues contribute disproportionately.
  • Interpretation: Target for optimization or mutation.

Usage Tips

  • Filter Settings:

    • Set min_contribution_percent to 5–10% to focus on strongest contributors.

  • Interpretation Scope:

    • Energies are not absolute binding free energies.
    • Use for relative comparison or pattern detection, not quantitative prediction.