RMSF Analysis: Reference Structure Selection

This guide explains how to choose the appropriate reference structure for RMSF (Root Mean Square Fluctuation) analysis in PolyzyMD. The choice of reference affects the interpretation of your results and should match your scientific question.

Overview

RMSF measures how much each residue fluctuates around a reference position during a simulation. Before computing RMSF, the trajectory must be aligned to a reference structure to remove global translation and rotation.

The choice of reference structure determines what you’re measuring:

Reference Mode	What RMSF Measures
`centroid`	Fluctuations around the equilibrium (most visited) conformation
`average`	Pure thermal fluctuations around the mathematical mean
`frame`	Fluctuations relative to a specific functional state
`external`	Deviations from a condition-independent external structure (e.g., crystal structure)

Reference Modes

Centroid (Default)

calc = RMSFCalculator(
    config,
    reference_mode="centroid",  # Default
    equilibration="100ns",
)

The centroid mode finds the most populated conformational state using K-Means clustering on all protein atoms. It returns the actual frame from the trajectory that best represents where the protein spends most of its time.

When to use:

Standard protein flexibility analysis
You want to know how much residues deviate from the equilibrium state
You have a well-behaved single-state protein

Implementation details:

Uses K-Means clustering with k=1 on all protein atoms (including side chains)
Finds the frame closest to the cluster center
This is a real sampled conformation, not a synthetic structure

Average Structure

calc = RMSFCalculator(
    config,
    reference_mode="average",
    equilibration="100ns",
)

The average mode computes the mathematical mean of all atomic positions and aligns the trajectory to this synthetic structure.

When to use:

You want a pure measure of thermal fluctuations
Comparing to literature values that used average-based RMSF
The protein samples a single conformational basin

Caveats:

The average structure may have unphysical geometry (distorted bond lengths, angles) because it’s a mathematical construct, not a real conformation
If the protein transitions between distinct states, the average may not represent any physically meaningful conformation

Specific Frame

calc = RMSFCalculator(
    config,
    reference_mode="frame",
    reference_frame=500,  # 1-indexed frame number
    equilibration="100ns",
)

The frame mode aligns to a user-specified frame. This is powerful for analyzing fluctuations relative to a known functional state.

When to use:

Enzyme catalysis studies where you want to measure deviations from a “catalytically competent” conformation
Comparing flexibility before/after a conformational change
You have identified a specific frame of interest (e.g., substrate bound, active site properly configured)

Example: Catalytic Competence Analysis

For enzyme studies, you might want to find a frame where:

The catalytic triad residues are properly hydrogen-bonded
The substrate is positioned correctly
The active site is in the “reactive” configuration

# First, identify the catalytically competent frame
# (e.g., by analyzing hydrogen bond distances, active site geometry)
catalytic_frame = 1523  # Frame where catalytic triad is aligned

# Then compute RMSF relative to that state
calc = RMSFCalculator(
    config,
    reference_mode="frame",
    reference_frame=catalytic_frame,
    equilibration="100ns",
)
result = calc.compute(replicate=1)

# High RMSF values now indicate residues that frequently deviate
# from the catalytically competent geometry

External PDB

Python

calc = RMSFCalculator(
    config,
    reference_mode="external",
    reference_file="/path/to/crystal_structure.pdb",
    equilibration="100ns",
)
result = calc.compute(replicate=1)

YAML (analysis.yaml)

rmsf:
  enabled: true
  selection: "protein and name CA"
  reference_mode: "external"
  reference_file: "/path/to/crystal_structure.pdb"

CLI

polyzymd analyze rmsf -c config.yaml -r 1 \
    --reference-mode external \
    --reference-file /path/to/crystal_structure.pdb \
    --eq-time 100ns

The external mode uses an external PDB file — typically a crystal structure — as both the alignment target and the reference for RMSF computation. Unlike the other modes, the reference is not derived from the trajectory itself, making it condition-independent.

RMSF then measures how much each residue deviates from the external structure’s geometry during the simulation:

\[ \text{RMSF}_i = \sqrt{\frac{1}{T} \sum_{t=1}^{T} \left( \mathbf{r}_i(t) - \mathbf{r}_i^{\text{ext}} \right)^2} \]

where \(\mathbf{r}_i^{\text{ext}}\) is the position of atom \(i\) in the external PDB, rather than the trajectory average \(\langle \mathbf{r}_i \rangle\).

When to use:

Enzyme catalysis studies where a crystal structure represents the catalytically competent geometry
Comparing how different conditions (polymer baths, temperatures, mutants) affect deviation from a known functional conformation
You need a condition-independent reference so RMSF values are directly comparable across conditions
Comparing to experimental B-factors derived from the same crystal structure

Validation: PolyzyMD validates that the external PDB’s protein atoms match the simulation topology exactly (same atom count for the selected atoms). If there is a mismatch, a clear error message reports the expected vs. actual atom counts.

Note

The external PDB must contain the same protein with the same residue numbering as your simulation. It does not need to contain solvent, ions, or polymer — only the protein atoms that match your selection string are compared.

Example: Catalytic Competence Analysis

For enzyme-polymer studies, a crystal structure provides a physically grounded reference for measuring how well each condition maintains the catalytically competent geometry:

# Use the crystal structure as the "gold standard" reference
calc = RMSFCalculator(
    config,
    reference_mode="external",
    reference_file="/path/to/enzyme_crystal.pdb",
    selection="protein and name CA and resid 5:175",  # Trim flexible termini
    equilibration="100ns",
)
result = calc.compute(replicate=1)

# Lower RMSF = closer to crystal structure = more catalytically competent
# Compare active site residues across conditions
for resid, rmsf in zip(result.residue_ids, result.rmsf_per_residue):
    if resid in [77, 133, 156]:  # Catalytic triad (e.g., LipA)
        print(f"Residue {resid}: {rmsf:.2f} Å from crystal")

Tip

Residue range truncation: Flexible N/C-terminal loops often dominate RMSF statistics and obscure active-site signals. Use a residue range selection (e.g., "protein and name CA and resid 5:175") to exclude terminal residues. No code changes are needed — MDAnalysis handles the selection natively.

Alignment Selection

Independently of the reference mode, you can control which atoms are used for the alignment superposition:

calc = RMSFCalculator(
    config,
    reference_mode="centroid",
    alignment_selection="protein and name CA",  # Default: backbone alignment
    centroid_selection="protein",  # Default: all protein atoms for clustering
)

Alignment selection (alignment_selection):

Controls which atoms are used to superimpose frames
Default is "protein and name CA" (alpha carbons only)
This removes global rotation/translation of the protein backbone

Centroid selection (centroid_selection):

Controls which atoms are used for K-Means clustering (centroid mode only)
Default is "protein" (all protein atoms)
Using all atoms captures side chain conformations in the clustering

Practical Recommendations

For General Flexibility Analysis

Use the default centroid mode:

polyzymd analyze rmsf -c config.yaml -r 1 --eq-time 100ns

This gives you flexibility relative to the equilibrium state, which is usually what you want for comparing flexibility across conditions or mutants.

For Enzyme Mechanism Studies

Consider using external mode with a crystal structure, or frame mode with a catalytically relevant frame from the trajectory:

# Option 1 (Recommended): External crystal structure as reference
# Best when you have a high-resolution crystal structure and want
# condition-independent comparison
calc = RMSFCalculator(
    config,
    reference_mode="external",
    reference_file="/path/to/crystal_structure.pdb",
    selection="protein and name CA and resid 5:175",
    equilibration="100ns",
)

# Option 2: Specific trajectory frame as reference
# Useful when the relevant conformation differs from the crystal structure
# (e.g., ligand-induced conformational change during simulation)
calc = RMSFCalculator(
    config,
    reference_mode="frame",
    reference_frame=catalytic_frame_id,
    equilibration="100ns",
)

RMSF values will then tell you which residues deviate from the active conformation, potentially disrupting catalysis.

For Comparing to Published Data

Check what reference the published study used:

Many older studies use average structures
Some use the first frame or crystal structure
If the study compares to crystallographic B-factors, use external mode with the same crystal structure PDB
Match the methodology for fair comparison

Technical Notes

Why Alignment Matters

Without alignment, RMSF measures the total displacement of atoms, which includes:

Global translation (protein drifting through the box)
Global rotation (protein tumbling)
Internal flexibility (what we actually want)

Alignment removes (1) and (2), leaving only the internal flexibility signal. This is why unaligned RMSF values can be ~50 Å while properly aligned values are typically 0.5-5 Å for a stable protein.

Memory Considerations

All alignment modes use in-memory trajectory alignment (in_memory=True). For very long trajectories, this may require significant RAM. Consider:

Using a larger equilibration time to reduce frames
Processing replicates sequentially rather than in parallel

Reproducibility

The result files store all alignment metadata:

reference_mode: Which method was used
reference_frame: The specific frame (1-indexed, or None for average)
reference_file: Path to the external PDB (when using external mode)
alignment_selection: Which atoms were used for superposition

This ensures your analysis is fully reproducible.