1. Introduction

Molecular dynamics (MD) simulation is a cornerstone of computational structural biology, enabling the study of protein conformational dynamics, folding, and ligand binding. Analysis of MD trajectories requires computing structural metrics such as RMSD (root mean square deviation), RMSF (root mean square fluctuation), contact maps, and principal component analysis (PCA).

Existing tools such as MDtraj (McGibbon et al., 2015) and MDAnalysis (Michaud-Agrawal et al., 2011) provide comprehensive functionality but require complex installation and have large dependency trees. For educational use, rapid prototyping, and deployment in constrained environments, a pure NumPy implementation is preferable.

We present MDAnalysisEngine, which implements all core MD analysis algorithms using only NumPy, SciPy, and matplotlib. The tool downloads real PDB structures from RCSB and generates synthetic Langevin dynamics trajectories for demonstration, making all results fully reproducible without any MD simulation software.

2. Methods

2.1 Structure Retrieval

PDB structures are downloaded from RCSB (https://files.rcsb.org/download/). Ca coordinates are parsed from ATOM records using pure Python string parsing.

2.2 Trajectory Generation

Synthetic trajectories are generated using a simplified Langevin dynamics model: at each step, Gaussian noise proportional to a B-factor profile is added, with a harmonic restoring force toward the equilibrium structure. Terminal residues (positions 1-5 and 72-76) have 3x higher B-factors than core residues (positions 21-50), reflecting known flexibility patterns.

2.3 RMSD (Kabsch Algorithm)

For each frame, we compute the optimal rotation matrix using the Kabsch algorithm (Kabsch, 1976): (1) center both structures; (2) compute the covariance matrix H = P^T Q; (3) SVD decomposition H = USV^T; (4) correct for reflection using det(V^T U^T); (5) compute RMSD after optimal superposition.

2.4 RMSF

Per-residue RMSF is computed from production frames (after equilibration): RMSF_i = sqrt(mean_t((r_i(t) - mean_r_i)^2)), where mean_r_i is the mean position over production frames.

2.5 Contact Map

Ca-Ca contact maps are computed using scipy.spatial.distance.cdist. A contact is defined as Ca-Ca distance < 8.0 A (configurable). Contact frequency is averaged over every 5th production frame for efficiency.

2.6 PCA

Trajectory PCA is performed by flattening each frame to a 3N-dimensional vector, centering, and computing SVD. The top 10 principal components are retained. Free energy is estimated as F = -kT ln(P), where P is the 2D histogram of PC1 vs PC2 coordinates.

2.7 H-Bond Proxies

Backbone H-bond occupancy is estimated from Ca-Ca distances: residue pairs (i, i+3) and (i, i+4) with Ca-Ca < 6.5 A in >50% of production frames are counted as persistent H-bond proxies.

3. Results

3.1 Ubiquitin (1UBQ) Analysis

Ubiquitin (76 residues) is a well-characterized protein with known flexibility at the C-terminus. Our analysis of a 500-frame trajectory (100 equilibration + 400 production) yields:

• Mean RMSD: 0.699 A (production), indicating a stable trajectory
• Most flexible residue: 75 (RMSF=1.711 A), consistent with C-terminal flexibility
• Most rigid residue: 34 (RMSF=0.226 A), in the hydrophobic core
• Persistent contacts: 321 (Ca-Ca < 8 A in >50% of frames)
• PC1: 14.9% of variance; PC1+PC2: 24.6%
• H-bond proxies: 37 persistent backbone H-bonds

3.2 Free Energy Landscape

The PC1-PC2 free energy landscape shows a single dominant conformational basin, consistent with ubiquitin's stable globular fold. Minor secondary basins correspond to C-terminal fluctuations.

4. Availability

MDAnalysisEngine is available at https://github.com/junior1p/MDAnalysisEngine under the MIT license. All results are reproducible by running python md_analysis_engine.py --pdb 1UBQ --n_frames 500.

References

1. Kabsch W (1976). A solution for the best rotation to relate two sets of vectors. Acta Crystallographica A.
2. McGibbon RT, et al. (2015). MDTraj: A Modern Open Library for the Analysis of Molecular Dynamics Trajectories. Biophysical Journal.
3. Michaud-Agrawal N, et al. (2011). MDAnalysis: A toolkit for the analysis of molecular dynamics simulations. Journal of Computational Chemistry.