apzubarev/kinetic-metric-protein-landscapes: ULTRAMETRIC ANALYSIS OF PROTEIN ENERGY LANDSCAPES BASED ON KINETIC METRICS AND MARKOV STATE MODELS

ULTRAMETRIC ANALYSIS OF PROTEIN ENERGY LANDSCAPES BASED ON KINETIC METRICS AND MARKOV STATE MODELS (MSM) Author: A. P. Zubarev Version: 1.12 (Fixed architecture: MD is launched only MAIN_NUM_SAMPLES times, files are saved with suffixes, ALPHA_ANALYSIS uses prepared trajectories. Table output is unconditional, added alpha indication for MAIN_ANALYSIS)

Method

DESCRIPTION AND PROGRAM

Purpose

This program implements a method for testing the hypothesis of ultrametric organization of protein energy landscapes. The method is based on analyzing transition kinetics between macrostates (attraction basins) using a kinetic metric constructed via spectral decomposition of the symmetrized Kramers transition rate matrix and the Mahalanobis distance. VERSION 1.12 FEATURES: Removed the redundant MAIN_ALPHA_SAMPLES parameter. The number of replicas for ALPHA_ANALYSIS is now exactly MAIN_NUM_SAMPLES. FUNDAMENTAL ARCHITECTURE FIX: Molecular dynamics (MD) and subsequent trajectory clustering are performed STRICTLY MAIN_NUM_SAMPLES times. Simulation 1 generates files: {base}_trajectory_1.dcd, {base}_1.pdb, {base}_energies_1.txt. Simulation 2 generates files: {base}_trajectory_2.dcd, {base}_2.pdb, {base}_energies_2.txt. ALPHA_ANALYSIS mode no longer starts MD. It takes already generated (or loaded from disk) MAIN_NUM_SAMPLES trajectories and applies all alpha values from ALPHA_RANGE to each of them. This is physically correct and saves hundreds of hours of computation. In FROM_FILES_ANALYSIS mode, the program automatically scans the folder for files with suffixes (_trajectory_1.dcd, _trajectory_2.dcd, etc.), determines the number of available samples, and performs analysis (MAIN and ALPHA) for each of them. FINAL TABLES ARE NOW ALWAYS PRINTED, even if VERBOSE = False. The program performs the following computational stages: MOLECULAR DYNAMICS (MD) STAGE Generation of an equilibrium conformational ensemble by Langevin dynamics using the OpenMM package. Supported implicit-solvent force fields are GB99dms and Amber14/OBC2. ROBUST STRUCTURE PREPARATION: Before MD, the structure is automatically cleaned: non-protein residues (ligands, ions, water) are removed, the longest polypeptide chain is selected, and missing atoms/hydrogens are added. This ensures compatibility with OpenMM force fields for arbitrary PDB entries. IN-MEMORY BUFFERING: Instead of writing directly to disk via DCDReporter, the program uses a custom MemoryReporter that accumulates coordinates and energies in Python lists. This solves the fundamental problem of "garbage" frames when NaNs occur: lists in RAM can be atomically truncated, whereas binary DCD files cannot be rolled back. A single DCD file is written to disk via mdtraj ONLY ONCE at the very end, after successful completion of the simulation. If the number of rollbacks exceeds the limit (MD_MAX_ROLLBACKS), the integration time step is automatically reduced and the simulation is restarted. MARKOV STATE MODEL (MSM) CONSTRUCTION STAGE The continuous conformational space is discretized by clustering trajectory frames into microstates using the k-centers algorithm in RMSD metric. The atom selection for RMSD calculation is determined by the parameter GEOMETRY_RMSD_MODE: 'ca': only C-alpha atoms (faster, reflects backbone folding). 'all': RMSD over all atoms (stricter, accounts for side chains). IMPORTANT: The selected metric is used CONSISTENTLY both for clustering and for constructing the neighborhood graph, ensuring internal consistency of space discretization. For each microstate (cluster), the free energy is computed through the partition function over energies of frames inside that cluster. The neighborhood graph is built by one of two methods: Energetic: two clusters are neighbors if the difference of their free energies does not exceed the threshold dF_max = -RT ln(eta). This criterion is invariant to energy shifts and provides correct connectivity at all energy levels. Adaptive topological graph: for each cluster, an effective radius is computed (the specified percentile of distances from the center to cluster points). Two clusters are neighbors if the distance between their centers does not exceed alpha * (R_i + R_j), where alpha is NEIGHBOR_ALPHA. Based on free energies and the neighborhood graph, the symmetrized transition rate matrix K_sym is built. The off-diagonal elements K_sym_ij are computed by the Kramers formula: exp(-|F_i - F_j| / (2RT)) for neighboring microstates, otherwise 0. Macrostates (attraction basins) are identified by gradient descent on the free energy landscape. Attractors (local minima and plateaus) are determined automatically; each microstate is assigned to the basin reached by gradient descent. BASIN_EQ_THRESHOLD sets the energy equivalence threshold below which states belong to the same plateau. IMPORTANT ALGORITHM NOTE: Attractors are identified first (states that have no neighbors strictly lower than F_i - epsilon_eq). The equivalence graph analysis is performed ONLY among attractors to correctly merge degenerate minima/plateaus. Gradient descent is applied exclusively to unlabeled slope states. The analysis is performed strictly by connected components. For each component with number of basins >= MIN_BASINS_PER_COMPONENT, ultrametricity is computed independently. The final ultrametricity is computed as a weighted average. The fragmentation index f_inter is also computed. SPECTRAL ANALYSIS AND ULTRAMETRICITY TESTING STAGE For each significant connected component, spectral decomposition of K_sym_c is performed. Noise modes are filtered by the spectral gap criterion. Basin characteristic vectors are projected onto physical eigenmodes. The distance between basins is the Mahalanobis distance: a weighted Euclidean norm of projection differences with weights 1/|lambda_k|. Slower relaxation modes contribute more to the distance. Ultrametricity is tested by enumerating all triplets of basins. A triplet is considered nontrivially ultrametric if the two largest distances are approximately equal (relative difference <= epsilon), and the third one is significantly smaller (relative difference > delta). The fraction of such triplets is u_nt. EXECUTION ARCHITECTURE AND MODES The program operates in one or more of the following modes: FROM_FILES_ANALYSIS: Skip MD, load existing DCD/CIF/PDB files. The program automatically finds all samples (_1, _2, ...), clusters them, and builds MSMs. Then the enabled analytical modes are run. SENSITIVITY_ANALYSIS: Run analysis for several BASIN_EQ_THRESHOLD values specified in SENSITIVITY_RANGE. Cannot be combined with ALPHA_ANALYSIS. ALPHA_ANALYSIS: Run analysis for several NEIGHBOR_ALPHA values specified in ALPHA_RANGE. Cannot be combined with SENSITIVITY_ANALYSIS. The number of replicas for each alpha is exactly equal to the number of loaded or generated trajectories (MAIN_NUM_SAMPLES). MAIN_ANALYSIS: Analysis with base parameters for all generated or loaded trajectories. Mean, standard deviation, and CV are computed to assess stability. TECHNICAL REQUIREMENTS Execution environment: Windows 10/11, Linux, macOS. Python 3.10-3.12. Required packages (pip install): openmm, pdbfixer, mdtraj, numpy, scipy, threadpoolctl, torch HARDWARE REQUIREMENTS An NVIDIA CUDA GPU is recommended for MD. Automatic fallback to CPU. MSM and spectral analysis stages are performed on CPU. USER PARAMETERS All configurable parameters are grouped below with detailed comments.