Scientific software · Molecular dynamics
MDSmooth
Abstract
A molecular-dynamics trajectory is reduced to a handful of key frames by low-pass filtering a single per-frame collective variable and keeping only the frames where it turns. A zero-phase Butterworth filter (filtfilt, forward and backward, so peak timing never shifts) strips thermal jitter from an RMSD, PC1, tICA IC1, or dihedral-PCA series; the local maxima and minima of the filtered curve become morph key frames, interpolated along ChimeraX's corkscrew screw path at per-segment step counts equal to the real frame spacing. The filter has one degree of freedom, the low-pass cutoff, solved by bisection to hit a target key-frame count. A kinetic mode instead clusters the slow tICA coordinates into metastable states with a Markov State Model and PCCA+ and tours one representative frame per state. Every signal reports a cosine content: above 0.85 the slow mode is undersampled diffusion masquerading as a transition, and the extrema the filter found are noise.
keywords: molecular dynamics · zero-phase filtering · tICA · Markov state model
## From trajectory to key frames
The object filtered is a per-frame scalar series, one value summarizing the structure at each frame. The default is best-fit RMSD to a reference frame over the macromolecule and any ligand, with bulk solvent and free ions dropped so a box of diffusing water does not swamp the signal. A zero-phase Butterworth low-pass filter of order 5 (scipy.signal.filtfilt, run forward then backward to cancel phase shift) removes the high-frequency jitter. The significant frames are the local maxima and minima of the filtered series (scipy.signal.find_peaks), plus the first and last: the turning points where the motion reverses. A new coordinate set interpolates between them, each segment taking a step count equal to its frame spacing so the morph keeps the simulation's true timing. Interpolation defaults to ChimeraX's corkscrew morph engine, run once per segment, which turns the rigid-body change between key frames into a screw motion instead of dragging atoms through straight lines. A dependency-free linear fallback exists and is used automatically when the morph engine is absent.
## One cutoff, solved for a frame count
The cutoff frequency fc in cycles per frame is the filter's only knob. The expected count of local extrema of the filtered series rises monotonically with fc (Rice's formula), so the cutoff that yields about a target number of key frames is found by bisecting the band (0, Nyquist). Each probe is one filter-and-count pass; the whole search is tens of milliseconds. The count is an integer step function of the cutoff, so the search returns the cutoff whose count is closest to the target and prefers the smaller, smoother cutoff on ties. The rationale is spectral: slow conformational motion sits in the low-frequency, high-power part of the spectrum and thermal jitter spreads into a flat high-frequency shelf, so a good cutoff sits in the valley between them. The cutoff can also be set directly, or as the frequency retaining a fraction of cumulative spectral power. On a trajectory of about 1600 frames, fc = 0.01 gives roughly 20 key frames and 0.024 gives roughly 50. The power-fraction route is numerically touchy: an RMSD series is dominated by its mean offset, so a fixed-fraction rule swings between no valid filter and keeping most of the noise. The frame target is the default because it re-fits itself to whatever atoms and signal it is handed.
## Slow collective variables
RMSD is direction-blind: distinct conformations can share a magnitude. Four collective variables track the direction of motion instead. PC1 is the leading principal component of the Kabsch-aligned trajectory, the single highest-variance collective motion, computed from the temporal (frame-by-frame) covariance so the cost scales with frame count rather than coordinate count. IC1 is the slowest time-lagged independent component, the coordinate whose autocorrelation survives longest across a lag τ, obtained from the generalized eigenproblem
where C(0) is the instantaneous covariance and C(τ) the lagged covariance of the aligned, mean-removed coordinates; the slowest mode is the eigenvector of largest λ. With more coordinates than frames C(0) is singular, and tICA then returns a different spurious slow mode at every lag, so the coordinates are first projected onto the leading PCA subspace holding 95% of the positional variance and a small ridge is added to C(0). Dihedral PCA runs the same decomposition on backbone phi/psi angles mapped to (cos, sin), which needs no structural alignment and handles the periodic wrap-around. DeepTICA is a neural slow variable for nonlinear modes; it runs in an isolated PyTorch/mlcolvar virtual environment as a subprocess, trains several seeds, and reports their agreement, because a learned CV is only trustworthy when independent seeds converge.
## Guarding against undersampling
A trajectory that never crosses its barrier of interest diffuses like a high-dimensional random walk, and a random walk projects onto its leading PCA/tICA mode as a near-perfect half-cosine. The cosine content of the mean-removed series measures exactly that resemblance:
which is 1 for a pure cos(πt/n) and near 0 for a signal with no slow drift. High c means the slow mode may be undersampled diffusion, not a sampled transition. The threshold is set at 0.85, deep in the diffusive regime, not at Hess's 0.5: a single genuine slow transition is also roughly cosine-shaped, and RMSD-to-a-reference naturally drifts to c near 0.8, so a low threshold cries wolf on good signals. It is a caution, never a gate. For IC1 the tool adds a second check: it recomputes the mode at τ/2, τ, and 2τ and correlates the shapes, because a real slow mode keeps its shape across lag and a spurious one flips around. The verdict is blunt: when the cosine content is high, the extrema the filter found are noise artifacts, and no amount of filtering recovers a transition the simulation never sampled.
## Kinetic mode: metastable states
Kinetic mode answers a different question: not "sample the main motion" but "show the distinct shapes the molecule occupies". It runs on tICA and ignores the signal keyword. The slow coordinates are clustered into many microstates with k-means; a Markov State Model is estimated on the largest kinetically-connected set; PCCA+ coarse-grains those into the requested metastable states; each state's representative is the frame with the highest membership. The tour visits states in nearest-neighbour order through the full tICA space rather than along IC1 alone, so states that overlap on IC1 but differ on other slow coordinates are not crushed together. It pays off only when the trajectory is long enough that the molecule visits each shape several times and transitions between them; otherwise the Markov State Model has no rates to estimate. The command stops with a message ("only N kinetically-connected microstates; not enough for an MSM") rather than returning arbitrary slices of noise.
Python, NumPy, SciPy (Butterworth filtfilt, find_peaks, generalized eigh), deeptime (MSM / PCCA+), PyTorch + mlcolvar (DeepTICA, isolated venv), ChimeraX 1.5 bundle API, Qt.
Liam Kozma · liam@liamkozma.com