α-zein (A8HNE1) — Molecular Dynamics Analysis Pipeline

A complete, modular, publication-quality analysis pipeline for a 100 ns GROMACS MD simulation of an α-zein protein (UniProt A8HNE1, 187 residues) built as an ab initio AlphaFold 3 model and simulated in an explicit cellular/solvent environment (CHARMM36, TIP3P, K⁺/Cl⁻).

Every analysis is a self-contained Python module that can be run independently or through run_all.py. Each writes CSV data, 300-dpi PNG and PDF figures, and feeds an automatically generated Markdown report.

---

1. Why this pipeline is pure-Python

The production system is large (~608k atoms) and the trajectory is ~4.4 GB, but the protein is only 2,912 atoms. The pipeline therefore:

1. Preprocesses once — streams the full trajectory a single time and writes a PBC-corrected, protein-only trajectory (data/protein.xtc, ~21 MB, ~200× smaller). All protein analyses run on this and are fast.
2. Needs no external binaries — no gmx, no standalone DSSP, no AmberTools:
   • energies are read straight from the .edr with panedr,
   • secondary structure uses mdtraj's built-in DSSP,
   • MM-GBSA scripts are generated (see note below).

---

2. Project layout

analysis/
├── config.py                 # single source of truth for paths & parameters
├── run_all.py                # orchestrator (runs everything, builds the report)
├── requirements.txt / environment.yml / README.md
│
├── mdkit/                    # shared utilities
│   ├── io_utils.py           #   universe loading, CLI, CSV, progress bars
│   ├── plotting.py           #   Nature-like style, save PNG+PDF @300 dpi
│   └── stats_utils.py        #   mean/median/std/95% CI, rolling, block avg
│
├── scripts/                  # one module per analysis (each has run(cfg)+CLI)
│   ├── preprocess.py         #   build protein-only trajectory (run first)
│   ├── rmsd.py  rmsf.py  rog.py  sasa.py  com.py  structural.py
│   ├── hbonds.py  salt_bridges.py  native_contacts.py  contact_map.py
│   ├── pca.py  dccm.py  clustering.py  secondary_structure.py
│   ├── ramachandran.py  rin.py  prolif_analysis.py  mmgbsa.py
│   ├── energies.py  convergence.py  statistics_summary.py
│   └── report.py             #   Markdown report generator
│
├── data/       # protein.xtc, protein.pdb, meta.json (preprocessing outputs)
├── results/    # all CSV outputs (+ cluster_structures/*.pdb)
├── figures/    # all figures (PNG + PDF, 300 dpi)
├── tables/     # summary/statistics tables, MM-GBSA scripts
└── report/     # analysis_report.md, mmgbsa_workflow.md

---

3. Installation

# option A: conda (recommended — reliable rdkit/mdtraj/MDAnalysis builds)
conda env create -f environment.yml
conda activate zein-md

# option B: pip (Python 3.11)
python -m pip install -r requirements.txt

---

4. Usage

Run the entire pipeline

python run_all.py \
    --top  step5_production.tpr \
    --traj step5_production.xtc \
    --edr  step5_production.edr \
    --output .

Omit the flags to use the defaults in config.py. Preprocessing runs automatically the first time; use --force-preprocess to rebuild it.

Run a single analysis

python scripts/rmsd.py            # backbone & Cα RMSD
python scripts/pca.py             # PCA + free-energy landscape
python scripts/secondary_structure.py

Useful run_all.py options

python run_all.py --list                 # list module keys
python run_all.py --only rmsd,rmsf,rog    # run a subset
python run_all.py --skip prolif_analysis  # skip the slow ProLIF step

---

5. Analyses included

Module	Outputs
rmsd	backbone & Cα RMSD + moving average, plateau
rmsf	per-residue RMSF, flexible-region highlighting
rog	radius of gyration, mean ± σ
structural	end-to-end, Dmax, shape anisotropy κ², convex-hull volume, elongation
sasa	total + per-residue SASA (Shrake–Rupley)
hbonds	protein–protein & protein–solvent H-bonds, occupancy, lifetime, donor/acceptor table
com	COM X/Y/Z, distance travelled, 3D trajectory
pca	covariance, eigenvalues, PC1/PC2, free-energy landscape
dccm	dynamic cross-correlation matrix
contact_map	contact occupancy + time-resolved difference map
clustering	RMSD clustering, optimal k (silhouette), medoid PDBs, populations
secondary_structure	DSSP timeline, helix/sheet/coil %
salt_bridges	salt-bridge count & occupancy
native_contacts	fraction of native contacts Q(t)
ramachandran	φ/ψ before vs after MD + ensemble density
rin	residue interaction network + centralities/hubs
prolif_analysis	protein–water hydration frequency + ProLIF fingerprint
mmgbsa	gmx_MMPBSA workflow generation + structural proxy
energies	potential/kinetic/total energy, T, P, density from .edr
convergence	RMSD plateau, block averaging, first-50 vs last-50 ns, RMSIP
statistics_summary	master stats table + Pearson/Spearman correlations
report	Markdown report with automated biological interpretation

---

6. Note on MM-GBSA

MM-GBSA/MM-PBSA computes a binding free energy between two partners. This is a single, ligand-free protein, so a binding ΔG is not physically defined. The mmgbsa module therefore generates a ready-to-run gmx_MMPBSA workflow (tables/mmgbsa_run.sh, tables/mmpbsa.in) for use once two groups are defined, and provides a force-field-independent per-residue structural contribution proxy instead. See report/mmgbsa_workflow.md.

---

7. Reproducibility

• All parameters (cutoffs, strides, window sizes, temperature) live in config.py.
• results/pipeline_summary.json records every module's numeric summary and per-module timings.
• Figures are written as both PNG and PDF at 300 dpi for direct manuscript use.