Machine Learning Revolution: Next-Gen Coarse-Grained MD Parameters for Small Protein Folding

Abstract

This comprehensive review explores the transformative impact of machine learning and advanced parameterization methods on coarse-grained molecular dynamics (CGMD) for small protein folding studies. We cover foundational principles of traditional CG approaches like Gō models and UNRES, detail cutting-edge ML-driven methods including deep learning force fields and Bayesian optimization, and provide practical guidance for troubleshooting parameterization challenges. Through comparative analysis of different CG strategies and their validation against all-atom simulations and experimental data, this article equips researchers and drug development professionals with the knowledge to implement accurate, computationally efficient protein folding simulations that bridge molecular mechanisms with biomedical applications.

Foundations of Coarse-Graining: From Basic Principles to Advanced Landscape Theory

Protein folding, the process by which a polypeptide chain self-assembles into its unique three-dimensional native structure, is one of the most fundamental processes in structural biology. While all-atom molecular dynamics (MD) simulations, which explicitly represent every atom in the protein and solvent, provide extremely high-resolution spatial and temporal data on folding processes, they face insurmountable challenges when applied to many biologically relevant proteins [1]. The core issue is timescale disparity: biological folding events can range from microseconds to minutes, whereas all-atom MD simulations are typically limited to microsecond timescales even with specialized supercomputers [2] [1]. This computational bottleneck is particularly problematic for larger, multidomain proteins which comprise the majority of proteins in eukaryotes, bacteria, and archaea [2].

Coarse-grained (CG) modeling addresses this challenge by reducing the number of degrees of freedom in the system, typically by representing multiple atoms with a single "bead" or interaction site. This simplification allows simulations to access dramatically longer timescales and larger systems while retaining the essential physics of the folding process. As the demand to study complex folding phenomena, misfolding diseases, and protein design grows, coarse-graining has transitioned from a specialized technique to an essential component of the computational biophysicist's toolkit [3].

The Computational Bottleneck: All-Atom Limitations

The limitations of all-atom molecular dynamics become apparent when examining the relationship between protein size, folding time, and computational feasibility. At a biological level, the median length of a protein is 532, 365, and 329 amino acids in eukaryotes, bacteria, and archaea, respectively [2]. These proteins fold on timescales ranging from microseconds to tens of minutes, creating a fundamental sampling problem for conventional simulation approaches.

Table 1: Representative Protein Folding Timescales and Computational Challenges

Protein	Length (residues)	Experimental Folding Time	All-Atom Simulation Challenge
Trpcage	20	~4 μs	Achievable with current resources [1]
Villin Headpiece	35	~4-7 μs	Borderline achievable [1]
Pin1 WW Domain	35	>50 μs	Very difficult [1]
α1-antitrypsin	394	Minutes	Currently impossible [2]

Technically, all-atom MD simulations must use femtosecond timesteps (10⁻¹⁵ seconds) to maintain numerical stability, meaning that simulating a single microsecond of real time requires one billion integration steps. Specialized hardware such as ANTON or massively distributed computing projects like Folding@home have pushed the boundaries of what's possible, but even these resources are typically limited to simulating proteins of approximately 100 amino acids with folding times up to milliseconds [2]. For larger proteins that fold through long-lived intermediates corresponding to deep local energy minima, all-atom simulations are likely to remain trapped in limited conformational spaces for the foreseeable future [2].

Coarse-Grained Methodologies: Theoretical Foundations

Coarse-grained models for protein folding can be broadly categorized into several classes based on their theoretical foundations and design principles. The common goal is to reduce the number of degrees of freedom while preserving the essential features of the protein's energy landscape.

Structure-Based Models

Structure-based models, particularly Gō models, operate on the principle of minimal frustration, which states that naturally occurring protein sequences have evolved to minimize energetic conflicts between interactions in the native conformation [2] [4]. These models encode the native structure of the protein directly into the potential energy function, primarily favoring native contacts while ignoring or weakly repelling non-native interactions [2]. The potential energy in these models typically decreases with increasing numbers of native contacts, making the number of native contacts an effective reaction coordinate for folding [2]. Gō models can be simulated using either Monte Carlo or molecular dynamics methods and are computationally inexpensive compared to all-atom simulations, making them particularly valuable for studying folding mechanisms and intermediate states [2] [4].

Knowledge-Based and Physics-Based Models

Beyond pure structure-based approaches, other coarse-graining strategies include knowledge-based methods that derive interaction parameters from statistical analysis of known protein structures, and physics-based methods that aim to capture fundamental physicochemical principles [3]. The MARTINI force field, for instance, is widely used for studying biomolecular interactions though it has limitations modeling intramolecular protein dynamics [5]. Other approaches like UNRES and AWSEM attempt to balance structural accuracy with transferability across different protein systems [5].

Machine-Learned Coarse-Grained Potentials

Recent advances incorporate machine learning to develop transferable coarse-grained force fields. These approaches use neural networks trained on large datasets of all-atom simulations to learn effective potentials that can accurately reproduce protein thermodynamics while accelerating dynamics by several orders of magnitude [5] [6]. For example, CGSchNet demonstrates how a single transferable potential can simulate multiple proteins with different topologies, predicting metastable states of folded, unfolded, and intermediate structures consistent with experimental observations [5].

Table 2: Comparison of Major Coarse-Grained Methodologies

Methodology	Resolution	Theoretical Basis	Strengths	Limitations
Gō Models	1-3 beads per residue	Native structure encoding	High computational efficiency; Excellent for folding mechanisms	Limited transferability; Requires known native structure
Elastic Network Models	1 bead per residue (Cα)	Harmonic approximation around native state	Very fast; Good for near-native dynamics	Poor for large conformational changes
MARTINI	~4 heavy atoms per bead	Thermodynamic parameterization	Good for membranes & mixtures; Transferable	Limited intramolecular dynamics
Machine-Learned (CGSchNet)	Configurable	Learned from all-atom MD data	Transferable; Good thermodynamics	Complex training; Computational cost

Experimental Protocols and Implementation

Structure-Based Gō Model Implementation

Principle: This protocol creates a simplified representation of a protein where interactions present in the native structure are favorable, while non-native interactions are repulsive or neutral. The method is particularly effective for studying folding mechanisms and pathways [2] [4].

Step-by-Step Protocol:

• Native Structure Preparation:

  • Obtain the atomic coordinates of the protein's native state from the Protein Data Bank (PDB).
  • Identify the residues forming the native contacts. A common approach is to consider two residues to be in contact if any of their non-hydrogen side-chain atoms are within a specified cutoff distance (typically 4.5-6.5 Å) [4].

• Coarse-Grained Mapping:

  • Represent the protein backbone using a single bead per amino acid, placed at the Cα atom position [4] [3].
  • Define the chain connectivity through harmonic bonds between consecutive Cα beads.

• Potential Energy Function Setup:

  • Bonded interactions: Apply harmonic potentials for virtual bonds and angles to maintain chain geometry. Dihedral angle potentials can be added to mimic backbone conformational preferences [4].
  • Native contacts: Define an attractive potential (commonly a Lennard-Jones-like potential or a simple well) for each identified native contact. The energy well depth is typically set to a uniform value (e.g., ε) for simplicity, though flavored models can use different energies based on residue type [4].
  • Non-native interactions: Implement a purely repulsive potential (e.g., soft-sphere or repulsive Lennard-Jones) for all other non-bonded pairs to prevent chain crossing and unrealistic compactness.

• Simulation Execution:

  • Use molecular dynamics or Monte Carlo methods to simulate the system.
  • Employ enhanced sampling techniques like replica exchange molecular dynamics (REMD) to improve sampling over energy barriers, especially for proteins with complex landscapes [2] [3].

• Analysis:

  • Use the fraction of native contacts (Q) and root-mean-square deviation (RMSD) as reaction coordinates to analyze folding pathways and identify intermediates [5].

Machine-Learned Coarse-Grained Potential

Principle: This protocol uses neural networks to learn a coarse-grained force field directly from all-atom molecular dynamics data, creating a transferable model that can simulate multiple proteins [5] [6].

Step-by-Step Protocol:

• Training Data Generation:

  • Perform extensive, unbiased all-atom MD simulations of a diverse set of proteins, capturing both folded and unfolded states. The dataset should include proteins with varied secondary structures (α-helical, β-sheet, mixed) [6].
  • Extract coordinate-force pairs from these trajectories.

• Coarse-Grained Mapping:

  • Define a mapping from all-atom to coarse-grained representation. A common choice is retaining only Cα atoms, but backbone or other mapped beads can also be used [6].

• Neural Network Potential Training:

  • Employ a force-matching loss function: Train the neural network (e.g., CGSchNet) to predict coarse-grained forces such that the sum of the learned force and prior forces matches the mapped all-atom forces [6]. The loss function is: L(θ) = (1/3nM) Σ∥ΞF(r_c) + ∇U(Ξr_c; θ)∥²
  • Incorporate prior knowledge: Add simple analytical potential terms (a "prior") to maintain chain connectivity (bonds), prevent steric clashes (repulsion), and enforce correct chirality (dihedrals). This stabilizes training and prevents unphysical configurations [6].

• Production Simulation:

  • Run Langevin dynamics or other appropriate integrators using the trained neural network potential to compute forces.
  • Use parallel tempering to ensure converged sampling of the equilibrium distribution for free energy calculations [5].

• Validation:

  • Compare the free energy landscapes, metastable states, and structural fluctuations (RMSF) of the CG simulations against reference all-atom data or experimental observations [5] [6].

Diagram Title: Machine-Learned CG Model Workflow

Table 3: Essential Resources for Coarse-Grained Protein Folding Studies

Resource Category	Specific Tool/Method	Function & Application
Structure-Based Models	Gō-like Models [2] [4]	Simulate folding mechanisms by favoring native contacts; ideal for pathway analysis.
Elastic Network Models	Dynamic ENM (DENM) [7]	Study near-native fluctuations and force-induced unfolding; uses Cα coordinates.
Transferable CG Force Fields	MARTINI [5] [3]	Model biomolecular interactions in membranes and complexes; not optimized for folding.
Transferable CG Force Fields	AWSEM [5]	Predict protein folding and conformational dynamics using knowledge-based potentials.
Machine-Learned CG	CGSchNet [5]	Transferable potential for predicting folded/unfolded states and mutant effects.
Enhanced Sampling	Replica Exchange MD (REMD) [3]	Accelerate conformational sampling by running parallel simulations at different temperatures.
Enhanced Sampling	Parallel Tempering [5]	Enhanced sampling technique to ensure convergence of equilibrium distributions.
Analysis Metric	Fraction of Native Contacts (Q) [5]	Reaction coordinate measuring progress toward native structure.
Analysis Metric	Cα Root-Mean-Square Deviation (RMSD) [5]	Measure of structural similarity to the native state.

Application Notes: From Theory to Practical Insights

The practical utility of coarse-grained models is demonstrated through numerous successful applications to complex folding problems. For large, topologically complex proteins like serpins (e.g., α1-antitrypsin), which have folding times as long as tens of minutes and a kinetically trapped functional state, coarse-grained simulations have provided critical insights into folding pathways and the structural characteristics of intermediates that are difficult to capture experimentally [2]. These simulations can inform hypotheses about how folding proteins interact with cellular quality control machinery [2].

Furthermore, coarse-grained models have proven valuable in studying the relationship between protein sequence and folding behavior. By introducing energetic heterogeneity into Gō models (e.g., by weighting contact energies according to physicochemical principles or statistical potentials), researchers have investigated how subtle sequence changes alter folding mechanisms. For example, such models have explained the different folding nucleation sites in proteins L and G, which share a common topology but differ in sequence, consistent with experimental findings [4].

Beyond single proteins, coarse-grained models are uniquely suited for studying protein folding and stability in biologically relevant crowded environments. Models that account for the intrinsic thermodynamics of folding and effective protein-protein interactions have revealed that macromolecular crowding can significantly impact protein stability, phase behavior, and aggregation propensity—phenomena critical to both cellular function and biopharmaceutical development [8].

Diagram Title: Coarse-Grained Methodology Selection Guide

Coarse-grained modeling has become an indispensable approach for tackling the multiscale challenge inherent in protein folding studies. By strategically reducing computational complexity, these methods enable researchers to access the biologically relevant timescales and system sizes that remain beyond the reach of all-atom molecular dynamics. From structure-based Gō models that illuminate folding mechanisms to emerging machine-learned potentials that offer unprecedented transferability, coarse-grained simulations provide a powerful framework for generating testable hypotheses about how proteins fold, misfold, and function. As these methodologies continue to evolve and integrate with experimental structural biology, they will play an increasingly vital role in advancing our understanding of protein dynamics and in the design of novel proteins for therapeutic and biotechnological applications.

Classical all-atom molecular dynamics (MD) simulations are a cornerstone of computational biology, providing atomistic detail on the structure and dynamics of biological macromolecules [9]. However, their capacity to capture biochemical processes is fundamentally limited by computational cost, restricting the accessible length and time scales far below those relevant for many cellular processes, such as protein folding or large complex assembly [9] [3]. To address this challenge, coarse-grained (CG) modeling emerged as a reductionist strategy, wherein groups of atoms are consolidated into single interaction sites, or "beads" [9]. This simplification drastically reduces the number of degrees of freedom, enables longer timesteps, and smooths the energy landscape, collectively accelerating dynamics by several orders of magnitude [9].

The development of CG models for proteins has a rich history, originating from the pioneering work of Levitt and Warshel [9] [3]. There is no unique way to coarse-grain a system, and the choice of methodology hinges on the desired level of granularity and the specific biological question at hand [9]. Early CG models varied widely in resolution, from one-bead-per-amino-acid representations to more detailed models incorporating backbone and side-chain chemical specificity [9]. This application note delineates the historical evolution, operational protocols, and inherent limitations of three foundational classes of traditional CG approaches: structure-based Gō models, the physics-based UNRES model, and the top-down Martini force field, framing the discussion within their application to small protein folding studies.

Traditional Coarse-Grained Approaches

Structure-Based Gō Models

Historical Context and Principles: Gō models are a class of structure-based models that utilize the known native structure of a protein to define its energy landscape [5]. They are grounded in the energy landscape theory, which posits that biomolecular folding funnels are minimally frustrated, meaning that interactions stabilizing the native state are dominant [5]. In these models, the energy landscape is biased toward a single, pre-defined native conformation, making them particularly suited for studying folding mechanisms and pathways for proteins whose native structure is known.

Methodology and Application: In a typical Cα-based Gō model, each amino acid is represented by a single bead located on its Cα atom [10]. The potential energy function is composed of local and non-local terms. The local terms maintain the chain's connectivity and stiffness, while the non-local terms are crucial for folding: attractive interactions are defined only for native contacts (amino acid pairs in close spatial proximity in the native state), and repulsive interactions are assigned to all other non-native pairs. This design ensures the native state is the global energy minimum. The AICG2+ model is a modern example of a structure-based CG model that employs such a potential energy function [10].

Key Limitations: The primary limitation of Gō models is their lack of chemical specificity and transferability. Since the energetic preferences are entirely derived from a single protein's native structure, these models cannot distinguish between the dynamic behaviors of different amino acid sequences that share a similar fold [9] [5]. They are inherently system-specific and cannot be used for predictive simulations of unknown protein structures, folding upon binding, or the study of sequences where non-native interactions play a significant role [5].

The UNRES Force Field

Historical Context and Principles: The UNited RESidue (UNRES) model is a more detailed, physics-based CG force field developed for protein folding studies [9] [5]. Unlike Gō models, UNRES is not structure-based; its parameters are not derived from a known native structure. Instead, it is designed to be a transferable force field that can, in principle, predict the folded structure of a protein from its sequence alone.

Methodology and Application: UNRES employs a higher resolution mapping scheme. Each amino acid is represented by two interaction sites: one for the peptide group and one for the side chain [9]. This allows UNRES to account for the chemical specificity of different side chains, which is critical for realistic packing interactions within and between proteins. The force field has been parameterized using experimental data and information from all-atom simulations, aiming to reproduce the correct thermodynamic behavior of proteins [9] [5].

Key Limitations: Despite its complexity and transferability aims, UNRES has shown limitations in consistently capturing alternative metastable states and complex folding pathways for certain proteins [5]. Its accuracy can be system-dependent, and it may not always correctly represent the balance between different secondary structure elements, such as the competition between helical and beta-sheet formation observed in some proteins [5].

The Martini Force Field

Historical Context and Principles: The Martini force field represents a distinct "top-down" philosophy in coarse-graining. It is a generic, widely applicable CG force field parameterized primarily against experimental thermodynamic data, such as oil/water partitioning coefficients and densities [9] [11] [12]. Originally developed for lipids, its success led to its extension to proteins, nucleic acids, and carbohydrates, making it one of the most popular CG models for biomolecular simulations [3] [12].

Methodology and Application: Martini typically maps on average four heavy atoms into a single CG bead, and several water molecules are represented by a single water bead [11]. The force field uses a standard functional form with Lennard-Jones and Coulomb potentials. Its key strength lies in accurately modeling intermolecular interactions, making it exceptionally good for studying protein-lipid interactions [12], membrane remodeling [11], and, as recently demonstrated, protein-ligand binding [12]. In binding studies, ligands are placed randomly in the solvent, and unbiased MD simulations can spontaneously sample binding events and pathways to buried pockets [12].

Key Limitations: A well-known shortcoming of the Martini model is its inaccurate modeling of intramolecular protein dynamics [5]. The force field is not well-suited for studying processes like protein folding, as it often fails to stabilize native protein structures from extended or unfolded conformations without the application of external restraints. Furthermore, the use of a CG water model can affect the representation of explicit solvation effects [11].

Table 1: Comparative Summary of Traditional Coarse-Grained Models

Feature	Gō Models	UNRES	Martini
Parameterization Philosophy	Structure-based (Bottom-up from a single native structure)	Physics-based (Bottom-up from AA data & experiment)	Thermodynamic-based (Top-down from experimental data)
Mapping Resolution	One bead per amino acid (Cα) [10]	Two beads per amino acid (peptide group & side-chain) [9]	~4 heavy atoms per bead [11]
Chemical Specificity	Low (based on native contacts, not chemistry)	High (considers side-chain identity) [9]	Moderate (limited set of bead types) [12]
Transferability	Not transferable (system-specific) [5]	Transferable (to new sequences) [9]	Highly transferable across biomolecules [12]
Primary Application	Protein folding pathways & mechanisms [5]	Protein structure prediction & folding [5]	Protein-lipid interactions, membrane systems, ligand binding [12]
Key Limitation	Cannot predict unknown structures or folding upon binding	Limited accuracy for alternative metastable states [5]	Poor description of intramolecular protein dynamics & folding [5]

Experimental Protocols

Protocol for Protein Folding with a Cα-Based Gō Model

This protocol outlines the setup and execution of a folding simulation for a small protein using a structure-based Cα Gō model, such as the AICG2+ model [10].

Research Reagent Solutions:

• Protein Structure: PDB file of the target protein's native state (e.g., 1FME for the BBA protein [13]).
• Simulation Software: A CG simulator that supports structure-based models (e.g., OpenCafeMol [10]).
• Contact Definition: A algorithm to define native atomic contacts from the PDB structure (e.g., using a heavy-atom distance cutoff of 4.5–6.0 Å).

Step-by-Step Methodology:

• System Setup:
  • Mapping: Extract the Cα coordinates from the native PDB file. Each Cα atom becomes the location of one CG bead representing an amino acid.
  • Topology Generation: Define the chain connectivity based on the protein sequence, creating bonds between consecutive Cα beads.
• Force Field Parameterization:
  • Bonded Terms: Define harmonic potentials for bonds and angles based on the native structure's geometry. Torsional potentials can be derived from statistical distributions or the native structure's dihedrals [10].
  • Non-Bonded Terms:
    • Identify all native contact pairs from the reference PDB structure.
    • For native contacts, assign an attractive potential (e.g., a 10-12 Lennard-Jones potential).
    • For all other non-native pairs, assign a repulsive potential (e.g., a soft repulsive potential or a standard 12-6 Lennard-Jones potential).
• Simulation Execution:
  • Energy Minimization: Minimize the energy of the initial structure (often the native state or an unfolded conformation) to remove any steric clashes.
  • Equilibration: Run a short simulation at the target temperature (e.g., 300K) to equilibrate the system.
  • Production Run: Perform a long molecular dynamics simulation (e.g., using Langevin dynamics for temperature control). Multiple replicates from different initial unfolded states are recommended to ensure adequate sampling of the folding landscape.
• Analysis:
  • Reaction Coordinates: Monitor the fraction of native contacts (Q) and the Cα root-mean-square deviation (RMSD) from the native structure.
  • Free Energy Landscape: Construct free energy surfaces as a function of Q and RMSD to identify metastable states, folding intermediates, and transition states.

Protocol for Protein-Ligand Binding with the Martini Force Field

This protocol describes how to set up an unbiased binding simulation to spontaneously sample ligand binding to a protein using the Martini force field, as demonstrated for T4 lysozyme and GPCRs [12].

Research Reagent Solutions:

• Protein and Ligand Coordinates: All-atom structures of the protein and ligand.
• CG Mapping Tools: Software for generating Martini CG topologies (e.g., martinize.py for proteins, cg.py for small molecules).
• Simulation Engine: A MD engine supporting the Martini force field (e.g., GROMACS).
• Solvation Box: A box of CG water beads (e.g., Martini water) and ions at physiological concentration.

Step-by-Step Methodology:

• Coarse-Graining:
  • Protein: Map the all-atom protein structure to its Martini CG representation. Apply elastic network model (ENM) restraints to maintain the protein's tertiary structure during the simulation, as intramolecular dynamics are not Martini's strength.
  • Ligand: Map the all-atom ligand structure to a CG representation. Parameterization is based on matching the ligand's partition free energy between water and octanol [12].
• System Assembly:
  • Place the CG protein in a simulation box large enough to accommodate ligand binding and unbinding events (e.g., a 10 nm cube for a small protein).
  • Solvate the system with CG water beads.
  • Add ions to neutralize the system and achieve the desired salt concentration.
  • Insert one or more ligand molecules randomly in the solvent, away from the protein.
• Simulation Execution:
  • Energy Minimization: Minimize the energy of the full system to remove bad contacts.
  • Equilibration: Equilibrate the system in the NVT and NPT ensembles to stabilize temperature and pressure.
  • Production Run: Conduct multiple long, unbiased MD simulations (e.g., 30 trajectories of 30 µs each for a total of 0.9 ms per ligand [12]). The use of multiple replicates increases the statistics for binding/unbinding events.
• Analysis:
  • Binding Pose: Calculate the RMSD of the bound ligand relative to an experimental reference structure (if available).
  • Ligand Density: Generate 3D density maps of the ligand around the protein to identify binding pockets and pathways.
  • Binding Free Energy: Compute the potential of mean force (PMF) along a reaction coordinate (e.g., distance from the binding site) and integrate to estimate the binding free energy.

Figure 1: A generalized workflow for setting up and running a coarse-grained molecular dynamics simulation, highlighting the initial critical choice of the CG model type.

The Scientist's Toolkit: Essential Materials and Reagents

Table 2: Key Research Reagents and Computational Tools for CG Simulations

Item Name	Function/Description	Application Context
CG Simulation Software (OpenCafeMol)	A specialized, GPU-accelerated MD simulator for residue-resolution CG models, offering significant speedups over CPU-based calculations [10].	Folding simulations with Gō-like models; membrane dynamics.
CG Simulation Software (GROMACS)	A versatile, high-performance MD engine that supports a wide range of force fields, including Martini.	All-purpose MD, commonly used for Martini simulations of protein-lipid systems and ligand binding.
Weighted Ensemble Sampling (WESTPA)	An enhanced sampling toolkit that runs multiple parallel replicas and resamples them based on progress coordinates to efficiently capture rare events [13].	Systematically exploring protein conformational landscapes and folding/unfolding transitions.
Elastic Network Model (ENM)	A method to apply harmonic restraints between CG beads that are close in space, used to maintain a protein's native fold during Martini simulations.	Stabilizing protein structure in Martini simulations focused on binding or assembly, not folding.
Martini Water Bead	A CG particle representing approximately 3-4 real water molecules, parameterized to reproduce bulk water properties and interfacial tension [11].	Solvation in Martini force field simulations.
Benchmark Protein Dataset	A curated set of small, fast-folding proteins (e.g., Chignolin, WW domain, BBA) with diverse topologies, used for validation [13].	Benchmarking and validating the accuracy of new or existing CG force fields for protein folding.

The historical development of coarse-grained models has provided researchers with a powerful set of tools to access biological scales that are unreachable by all-atom simulations. The traditional approaches—Gō models, UNRES, and Martini—each embody a different philosophical approach to coarse-graining, leading to distinct strengths and application niches. Gō models excel at elucidating folding pathways for known structures, UNRES aims for predictive folding from sequence, and Martini dominates the study of biomolecular interactions in complex environments like membranes.

However, their limitations are significant. The lack of transferability in Gō models, the occasional inaccuracy of UNRES for metastable states, and the inability of Martini to model protein folding highlight a longstanding challenge: the difficulty in creating a universal, transferable CG model that retains chemical specificity and can accurately predict protein thermodynamics and kinetics across a wide range of sequences and environments [9] [5]. A primary hindrance has been the efficient representation of multi-body interactions essential for realistic implicit solvation and protein thermodynamics [5].

The future of the field, as evidenced by recent literature, is being shaped by the integration of machine learning (ML). New ML-based CG force fields, such as CGSchNet, are now being developed using bottom-up approaches trained on diverse all-atom simulation data [5] [14]. These models leverage deep learning to capture complex, many-body interactions, demonstrating promising results in achieving chemical transferability, predicting metastable states, and calculating relative folding free energies for unseen protein sequences, all while maintaining a computational efficiency that is orders of magnitude greater than all-atom MD [5]. The emergence of such models, along with standardized benchmarking frameworks [13], suggests that the long-sought goal of a universal and predictive coarse-grained model for protein simulation may soon be within reach.

The study of protein folding, the process by which a polypeptide chain attains its functional three-dimensional structure, is a fundamental challenge in molecular biology. Energy landscape theory provides a powerful conceptual framework for understanding this process, portraying it as a navigation across a multidimensional surface where the topographical features—energy basins, barriers, and pathways—dictate folding kinetics and thermodynamics [15]. For small protein folding studies, where atomistic simulations often struggle to access biologically relevant timescales, coarse-grained (CG) modeling offers a practical computational approach by reducing the number of degrees of freedom in the system [16]. The effectiveness of any CG model hinges upon how accurately it represents the underlying energy landscape of the protein, making energy landscape theory the indispensable foundation for developing and parameterizing these simplified representations.

The fundamental connection between these concepts is established through the potential of mean force (PMF), which serves as the rigorous physical basis for coarse-grained force fields [16]. The PMF is the free energy surface obtained by statistically averaging over the fine-grained degrees of freedom that are omitted in the CG model. Consequently, motion on the CG landscape is governed not only by the mean forces derived from the PMF but also by friction and stochastic forces arising from the integrated secondary degrees of freedom, making Langevin dynamics a natural choice for describing CG molecular dynamics [16]. This theoretical foundation ensures that CG models can extend simulation capabilities to large time and size scales while preserving the essential thermodynamics of the folding process [6].

Theoretical Foundations

The Potential of Mean Force as the Physical Basis for CG Models

In coarse-grained molecular dynamics, the motion of CG sites is governed by the equation derived from statistical mechanical principles:

where the first term represents the mean force acting on the ith coarse-grained site, derived from the potential of mean force W(R); the second term corresponds to friction forces that depend on the time correlation of force fluctuations; and the final term represents stochastic forces [16]. The potential of mean force W(R) is defined as:

where V(R;r) is the potential energy function of the full system, R denotes the CG degrees of freedom, and r represents the fine-grained degrees of freedom that are averaged out [16]. This formulation ensures that the equilibrium distribution of the CG system matches that of the all-atom system in principle, providing a thermodynamically consistent foundation.

In practical implementations, the complex memory effects of the friction term are often simplified using Markovian assumptions, leading to the more familiar Langevin equation:

where G is the inertia matrix, U is the effective CG energy function, Γ is the friction tensor, and frand represents random forces [16]. This equation forms the basis for most practical CG-MD implementations in protein folding studies.

Funneled Landscapes and Biomolecular Foldability

Energy landscape theory introduces the crucial concept of the folding funnel—an energy landscape that is biased toward the native structure [17]. A well-designed CG model for protein folding should reproduce this funneled characteristic, where the energy decreases as the protein approaches its native conformation. This organization minimizes frustration in the folding process and explains how proteins can rapidly find their unique native states among astronomically many possible conformations. For CG models, this implies that the effective energy function must be designed to stabilize the native structure while providing minimal frustration in the landscape, a principle that has been successfully implemented in structure-based models [18] and more recently in machine-learned potentials [5].

Methodological Approaches for CG Model Development

Traditional Parameterization Methods

Traditional approaches to developing CG force fields for protein folding have employed various strategies to capture the essential features of the energy landscape. Structure-based models (or Gō-models) use the known native structure to define attractive interactions for native contacts, effectively creating a perfectly funneled landscape [18]. These models have been particularly successful in elucidating folding mechanisms but require a priori knowledge of the native state.

For more general applications, physics-based parameterization incorporates non-native interactions through empirical potentials. A typical CG energy function for protein-protein association studies includes:

where E1 and E2 are structure-based potentials for individual proteins, and E12 contains non-native electrostatic and hydrophobic interactions between proteins [18]. The electrostatic component typically uses a Coulomb potential with distance-dependent dielectric screening, while hydrophobic interactions are often modeled with Lennard-Jones-type potentials parameterized using statistical energies from databases like Miyazawa-Jernigan [18].

Machine Learning Approaches

Recent advances have introduced machine learning methods for developing thermodynamically consistent CG potentials. These approaches use neural networks trained on all-atom simulation data to capture the many-body interactions essential for realistic protein thermodynamics [6] [5]. The variational force-matching method trains neural network potentials (NNPs) by minimizing the loss function:

where Ξ is the mapping operator from all-atom to CG coordinates, F(rc) are all-atom forces, and U is the CG potential parameterized by θ [6]. This approach ensures that the CG forces match the projected all-atom forces, leading to models that can preserve the thermodynamics of the original system while accelerating dynamics by several orders of magnitude [6].

To enhance transferability across different protein sequences, modern ML-CG models incorporate prior physical knowledge through bonded and repulsive terms to maintain chain integrity, and chirality constraints to prevent mirror-image structures [6]. The resulting models demonstrate remarkable capability in simulating proteins not included in the training set, with low sequence similarity (16-40%), indicating that they learn general physical principles of protein folding rather than merely memorizing training examples [5].

Table 1: Comparison of CG Model Parameterization Approaches

Method	Theoretical Basis	Required Input	Transferability	Key Applications
Structure-Based Models	Energy landscape theory, funneled landscapes	Native structure	System-specific	Folding mechanisms, pathways [18]
Physics-Based Parameterization	Potential of mean force, statistical potentials	Amino acid sequence, physicochemical data	Moderate to high	Protein-protein association, folding [18]
Machine Learning CG	Variational force-matching, neural networks	All-atom MD trajectories	High (sequence transferable)	Multiple proteins with single potential, mutant predictions [6] [5]
Self-Learning Multiscale	Iterative AA/CG optimization	AA force field, initial CG model	System-specific	De novo structure prediction [19]

Application Notes: Case Studies in Small Protein Folding

Coarse-Grained Simulations of Fast-Folding Proteins

Recent machine learning approaches have demonstrated remarkable success in simulating the folding of small, fast-folding proteins. In one notable study, researchers created a unique dataset of unbiased all-atom molecular dynamics simulations totaling approximately 9 milliseconds for twelve different proteins with diverse secondary structures [6]. This dataset was used to train neural network-based CG potentials that could accelerate dynamics by more than three orders of magnitude while preserving the thermodynamics of the systems.

The trained CG models successfully identified relevant structural states in the ensemble with comparable energetics to all-atom systems. Importantly, the study showed that a single CG potential could integrate all twelve proteins and even capture experimental structural features of mutated proteins not included in the training, demonstrating significant transferability [6]. This approach represents a major advancement toward universal CG models that can be applied across protein families without reparameterization.

Transferable Model Performance on Diverse Targets

Further development of transferable CG models has shown excellent performance across a range of small protein folding applications. These models successfully predict metastable states of folded, unfolded, and intermediate structures, fluctuations of intrinsically disordered proteins, and relative folding free energies of protein mutants [5]. Key demonstrations include:

• Folding of Chignolin and TRP-cage: The CG model correctly stabilized the native states of these small fast-folding proteins, including the same misfolded state with misaligned tyrosine residues as found in reference atomistic simulations for chignolin [5].
• Beta-Beta-Alpha (BBA) protein: While the model captured folded and unfolded states, the relative free energy differences showed some discrepancies from reference simulations, highlighting the challenge of mixed α/β motifs [5].
• Larger protein systems: The approach successfully folded the 54-residue engrailed homeodomain (1ENH) and 73-residue de novo designed protein alpha3D (2A3D), with Cα root-mean-square fluctuations in the folded state closely matching all-atom references [5].

Table 2: Performance of ML-CG Models on Representative Protein Folding Tasks

Protein System	Size (residues)	Key Results	Comparison to Reference
Chignolin	10	Correct native state stabilization, misfolded state identification	Matches all-atom MD metastable states [5]
TRP-cage	20	Folding to native state with correct topology	Native state as global free energy minimum [5]
Villin Headpiece	35	Robust folding to native structure	Comparable folded state stability [5]
BBA Protein	28	Captured folded state as local minimum	Some discrepancy in relative state energies [5]
Engrailed Homeodomain	54	Correct folding from extended state, terminal flexibility	Similar fluctuation profile to AA simulations [5]

Experimental Protocols

Protocol: Machine Learning CG Potential Development

This protocol outlines the key steps for developing a transferable coarse-grained potential for protein folding studies using machine learning approaches based on recent successful implementations [6] [5].

Data Preparation and All-Atom Reference Simulations

• Protein Selection: Curate a diverse set of small proteins (typically 10-80 residues) with varied secondary structure elements (α-helices, β-sheets, mixed α/β).
• All-Atom MD Simulations: Perform extensive unbiased all-atom molecular dynamics simulations starting from disordered configurations using explicit solvent models. For each protein, aim for sufficient sampling to observe multiple folding/unfolding events (typically microsecond to millisecond timescales per protein).
• Trajectory Processing: Extract coordinates and forces from the all-atom trajectories at regular intervals (e.g., every 100 ps). The dataset should comprise millions of coordinate-force pairs across all simulated proteins.

Coarse-Grained Mapping and Prior Potential Definition

• Mapping Selection: Define the CG representation, typically one to two beads per amino acid residue. Common choices include Cα-only or Cα/Cβ representations.
• Prior Potential Specification: Implement analytical prior potentials to maintain chain geometry and prevent unphysical configurations:
  • Bond potential: Harmonic restraints between adjacent beads: E_bond = ∑k_b(r - r_0)²
  • Angle potential: Harmonic restraints for chain angles: E_angle = ∑k_θ(θ - θ_0)²
  • Repulsive potential: Soft core repulsion to prevent bead overlap: E_rep = ∑ε(σ/r)¹²
  • Chirality potential: Dihedral restraints to enforce correct chirality: E_dih = ∑k_φ[1 + cos(n(φ - φ_0))]

Neural Network Potential Training

• Architecture Selection: Choose a neural network architecture suitable for molecular potentials (e.g., SchNet, Behler-Parrinello networks).
• Training Setup: Implement the variational force-matching loss function: L(θ) = (1/3nM)∑||ΞF_AA(r_c) + ∇U_NN(Ξr_c;θ)||² where Ξ projects all-atom coordinates to CG space, FAA are all-atom forces, and UNN is the neural network potential.
• Model Training: Train the network using stochastic gradient descent or Adam optimizer with appropriate batch sizes and learning rate scheduling until convergence of the loss function.

Model Validation and Testing

• Within-Test Evaluation: Assess performance on proteins from the training set by comparing free energy surfaces, folded state stability, and dynamics to all-atom references.
• Transferability Testing: Evaluate on completely unseen proteins with low sequence similarity to training examples (<40% identity).
• Experimental Comparison: Where possible, compare simulation results with experimental data such as folding free energies, mutant effects, and structural fluctuations.

Protocol: Enhanced Sampling for Protein Folding Landscapes

This protocol describes the push-pull-release (PPR) sampling method for efficient exploration of protein folding and association landscapes, adapted from successful implementations in CG protein studies [18].

System Setup and Initialization

• CG Model Preparation: Initialize the CG representation of the protein system using either a physics-based or machine-learned potential as described in Section 5.1.
• Simulation Parameters: Set up Langevin dynamics simulations at physiological temperature (300 K) with an appropriate friction coefficient (typically 50 ps⁻¹ for CG simulations) and time step (0.01 ps for CG-MD).
• Initial Configuration Generation: For association studies, generate multiple initial configurations by random translation and rotation of protein partners to ensure diverse starting orientations.

Push-Pull-Release Cycling

• Pull Phase: When proteins are in close association (e.g., minimum residue-residue distance < 7.6 Å), apply a biasing potential to pull them apart: E_PPR = k(R - R_t)² where R is the instantaneous center-to-center distance and R_t is a target separation distance.
• Push Phase: When proteins are separated beyond a threshold distance (R_c), apply a biasing potential to push them together using the same functional form.
• Release Phase: When proteins are within interaction distance but not yet in a deep energy minimum, remove the biasing potential and allow unbiased interaction for a predefined simulation period (e.g., 10 ns).
• Cycle Repetition: Repeat the push-pull-release cycle multiple times (typically hundreds to thousands of cycles) to thoroughly sample the association and folding landscape.

Trajectory Analysis and State Identification

• Configuration Clustering: Apply two-step clustering to identify metastable states:
  • Orientation clustering: Align structures based on one protein and cluster using key Cartesian coordinates from the partner.
  • RMSD clustering: Further cluster configurations using pairwise RMSD-based clustering with an appropriate cutoff (e.g., 5 Å).
• Energy Landscape Projection: Project simulation data onto relevant collective variables (e.g., fraction of native contacts, RMSD, radius of gyration) to construct free energy surfaces.
• State Characterization: For each identified metastable state, analyze structural features, interaction patterns, and relative populations.

The Scientist's Toolkit

Table 3: Essential Research Reagents and Computational Tools for CG Protein Folding Studies

Tool/Resource	Type	Function	Example Applications
All-Atom Force Fields	Software parameter set	Provides reference data for CG model development	CHARMM, AMBER, OPLS for reference MD [6]
CG Simulation Packages	Software	Performs coarse-grained molecular dynamics	GROMACS, LAMMPS, CHARMM with CG capabilities [18]
Variational Force-Matching	Algorithm	Trains ML-CG potentials to match all-atom forces	Developing transferable protein potentials [6]
Structure-Based Models	CG potential	Creates funneled landscapes using native structure	Folding mechanism studies [18]
Push-Pull-Release Sampling	Enhanced sampling method	Accelerates exploration of binding/unbinding events	Protein-protein association studies [18]
Neural Network Potentials	ML architecture	Captures multi-body interactions in CG models	Transferable CG force fields [5]

Energy landscape theory provides the essential conceptual foundation for developing effective coarse-grained models of protein folding. By representing the folding process as navigation on a multidimensional surface, this theory guides the design of CG models that capture the essential thermodynamics and kinetics while offering computational efficiency orders of magnitude greater than all-atom approaches. Recent advances in machine learning, particularly through variational force-matching and neural network potentials, have enabled the development of transferable CG models that maintain physical accuracy while extending simulation capabilities to biologically relevant timescales and system sizes. The continued integration of energy landscape principles with data-driven approaches promises to further advance our ability to simulate and understand protein folding, with significant implications for protein design and therapeutic development.

Coarse-grained (CG) molecular dynamics simulations have emerged as an indispensable tool in the study of biomolecular systems, particularly for investigating large and complex processes like protein folding that remain challenging for all-atom methods. By representing groups of atoms as single interaction sites or "beads," CG models significantly reduce the computational cost of simulations, enabling researchers to study larger systems for longer timescales [20]. This approach is firmly grounded in statistical mechanical descriptions of protein dynamics, such as energy landscape theory, which suggest that a protein's free energy landscape can be sufficiently described by a reduced number of collective variables with minimal loss of accuracy compared with atomistic molecular dynamics [5].

The development of accurate and transferable CG force fields represents an ongoing challenge in computational biophysics. The primary hindrance has been the difficulty in efficiently modeling multi-body interaction terms, which are essential to realistically represent correct protein thermodynamics and implicit solvation effects [5]. However, recent advances in deep learning have enabled the machine-learning of such many-body CG force fields, bringing the field closer to universal, quantitative, and predictive CG models that rival the reliability of modern atomistic force fields [5].

This application note focuses on three fundamental CG representations: Cα models, multi-bead schemes, and implicit solvent approaches, providing researchers with detailed protocols and comparisons to facilitate their implementation in small protein folding studies.

Key Coarse-Grained Representations

Cα Models

Cα models represent one of the most simplified yet powerful approaches to coarse-graining proteins. In these models, each amino acid residue is represented by a single bead, typically located at the Cα carbon position. This drastic reduction in degrees of freedom allows for unprecedented computational efficiency in simulating protein folding and dynamics.

The key advantage of Cα models lies in their ability to capture the essential structural and dynamic properties of proteins while being several orders of magnitude faster than all-atom models [5]. Recent machine-learned Cα models have demonstrated remarkable capabilities, including predicting metastable states of folded, unfolded, and intermediate structures, fluctuations of intrinsically disordered proteins, and relative folding free energies of protein mutants [5].

Table 1: Characteristics of Cα Coarse-Grained Models

Feature	Description	Common Applications
Resolution	One bead per amino acid, typically at Cα position	Protein folding, large-scale conformational changes
Computational Speed	Several orders of magnitude faster than all-atom MD [5]	Long-timescale simulations, enhanced sampling
Accuracy	Predicts folded states with low Cα RMSD and high fraction of native contacts [5]	Folding mechanisms, metastable state identification
Transferability	Machine-learned versions show promise for transferability across sequences [5]	Simulation of proteins not in training set
Common Force Fields	Structure-based Gō model, AICG2+ [21], machine-learned potentials [5]	Fast-folding proteins, intrinsically disordered proteins

Multi-Bead Schemes

Multi-bead schemes offer a more detailed representation than Cα models by incorporating additional beads to represent key structural elements of amino acids. These models typically include beads for the backbone and various sidechain components, providing a more nuanced description of molecular interactions.

Popular multi-bead frameworks include the Martini force field, which typically maps 4-5 heavy atoms to a single CG bead, and SIRAH, which employs a hybrid resolution approach with different levels of coarse-graining for protein backbone and side chains [20] [22]. The MARTINI force field has been particularly widely adopted, though recent developments have addressed its limitations in modeling intramolecular protein dynamics [5].

These models excel in capturing specific interactions that drive biological processes. For instance, the HPS model for intrinsically disordered proteins (IDPs) and the 3SPN series models for nucleic acids have provided valuable insights into biomolecular condensation and genome organization [21]. The multi-bead representation allows for a more accurate description of chemical specificity compared to single-bead approaches, while maintaining significant computational advantages over all-atom models.

Implicit Solvent Approaches

Implicit solvent models dramatically enhance computational efficiency by eliminating explicit solvent degrees of freedom, instead representing solvent effects through effective interactions and mean-field approximations. These approaches can be derived through rigorous statistical mechanical frameworks, such as the Mori-Zwanzig formalism, which links CG variables directly to microscopic dynamics [23].

In these models, the dynamics of the CG system is governed by the generalized Langevin equation (GLE), where solvent-mediated effects are captured through non-Markovian stochastic dynamics [23]. The memory kernel, which can be computed from all-atom molecular dynamics trajectories, encapsulates the friction and random forces exerted by the omitted solvent molecules. To circumvent the difficulty in direct evaluation of the memory term and generation of colored noise, the non-Markovian dynamics can be mapped to Markovian dynamics in an extended space by supplementing the CG system with auxiliary variables [23].

Table 2: Comparison of Implicit Solvent Methodologies

Method	Theoretical Basis	Implementation	Advantages
Generalized Langevin Equation	Mori-Zwanzig formalism [23]	Auxiliary variables coupled to momentum	Conserves dynamic properties of reference system
Effective Energy Terms	Mean-field approximations	Modified non-bonded interactions	Computational efficiency, simplified analysis
Machine-Learned Solvation	Deep learning on all-atom simulation data [5]	Neural network potentials	Transferability, many-body effects

Quantitative Comparison of CG Representations

The selection of an appropriate CG representation depends on the specific research question, system size, and properties of interest. The table below provides a comprehensive comparison of the key characteristics across different CG approaches.

Table 3: Quantitative Comparison of Coarse-Grained Representations

Parameter	Cα Models	Multi-Bead Schemes	Implicit Solvent
Beads per Residue	1	2-5 (backbone + sidechains)	1-5 (depending on base model)
Speedup vs All-Atom	Several orders of magnitude [5]	10-1000x	Additional 10-100x beyond explicit solvent CG
Accuracy (RMSD)	~0.5 nm for 54-residue engrailed homeodomain [5]	Higher than Cα models	Comparable to underlying structural model
Memory Requirements	Lowest	Moderate	Low to moderate
Best Applications	Large-scale dynamics, folding pathways	Protein-protein interactions, membrane systems	Large systems, slow processes
Limitations	Limited chemical specificity	Parameterization complexity	May miss specific solvent effects
Example Force Fields	AICG2+ [21], machine-learned [5]	MARTINI [20] [22], SIRAH [22], ELBA [20]	GLE with auxiliary variables [23]

Research Reagent Solutions

Successful implementation of CG-MD simulations requires a suite of specialized software tools and force fields. The table below catalogues essential "research reagents" for conducting state-of-the-art CG simulations.

Table 4: Essential Research Reagents for Coarse-Grained Simulations

Reagent	Type	Function	Application Context
GENESIS CGDYN [21]	Software	MD engine with dynamic load balancing	Large, heterogeneous biomolecular systems
GROMACS [20] [22]	Software	High-performance MD package	General biomolecular simulations
MARTINI 3 [22]	Force Field	Multi-bead CG force field	Biomolecular interactions, membranes
SIRAH [22]	Force Field	Multi-bead CG force field	Proteins, DNA, and their complexes
AICG2+ [21]	Force Field	Structure-based model for proteins	Protein folding, dynamics
HPS Model [21]	Force Field	IDP-specific potential	Intrinsically disordered proteins, LLPS
3SPN Series [21]	Force Field	Nucleic acid models	DNA, RNA structure and dynamics
Machine-Learned Potentials [5]	Force Field	Transferable bottom-up potentials	Multiple proteins with sequence transferability

Experimental Protocols

Protocol for Machine-Learned Cα Model Simulation

This protocol outlines the steps for implementing machine-learned Cα model simulations based on the CGSchNet framework [5], which has demonstrated success in predicting protein folding landscapes and metastable states.

• Training Data Generation: Generate a diverse dataset of all-atom explicit solvent simulations of small proteins with varied folded structures, including combinations of dimers of mono- and dipeptides [5].
• Model Training: Train the CG force field (CGSchNet) using the variational force-matching approach on the reference all-atom simulation data to learn effective interaction potentials [5].
• System Setup: For new protein sequences, initialize structures from extended configurations or known experimental structures.
• Equilibration: Perform structural relaxation simulations to eliminate bad contacts and reach stable thermodynamic states.
• Production Simulation: Conduct parallel-tempering (PT) simulations or long constant-temperature (300 K) Langevin simulations to ensure converged sampling of the equilibrium distribution [5].
• Validation: Compare key metrics such as fraction of native contacts (Q) and Cα root-mean-square deviation (r.m.s.d.) against reference all-atom simulations or experimental data [5].
• Analysis: Calculate free energy surfaces using collective variables such as RMSD and fraction of native contacts to identify metastable states and folding pathways.

Figure 1: Cα Model Simulation Workflow

Protocol for Implicit Solvent CG-MD with GROMACS

This protocol provides a framework for conducting implicit solvent CG-MD simulations using GROMACS, adaptable for studying both folded and disordered proteins and peptides [20] [22].

• System Preparation: Obtain protein structure from PDB or generate using modeling software. For the amyloid-β42 example, structures would be prepared for both wild-type and disulfide-bonded variants [22].
• Conversion to CG Representation: Convert the all-atom structure to a coarse-grained model using appropriate mapping. For MARTINI, this typically involves mapping 4-5 heavy atoms to a single CG bead [20].
• Force Field Selection: Choose an appropriate CG force field (MARTINI 3, SIRAH, etc.). For IDPs, consider scaling protein-water interactions to avoid over-compactness [22].
• Simulation Box Setup: Place the CG protein in an appropriate simulation box with sufficient padding for the intended dynamics.
• Energy Minimization: Perform energy minimization to remove steric clashes and bad contacts using steepest descent or conjugate gradient algorithms.
• Equilibration: Conduct equilibration simulations with position restraints on protein heavy atoms to relax the solvent and ions around the protein.
• Production Run: Execute production MD simulation with an integration time step appropriate for the CG model (typically 10-30 fs). For enhanced sampling, consider replica exchange MD [22].
• Trajectory Analysis: Analyze trajectories using GROMACS tools or custom scripts to calculate properties such as radius of gyration, RMSD, RMSF, and secondary structure content [20].

Figure 2: Implicit Solvent CG-MD Workflow

Advanced Applications and Case Studies

Protein Folding Landscape Mapping

Machine-learned Cα models have demonstrated remarkable capability in mapping protein folding landscapes. Recent work has shown that these models can successfully predict metastable states of folded, unfolded, and intermediate structures for various fast-folding proteins, including chignolin, TRPcage, BBA, and the villin headpiece [5]. The free energy landscapes generated from these simulations closely match reference all-atom simulations, with correctly populated native states characterized by high fraction of native contacts (Q close to 1) and low Cα root-mean-square deviation values [5].

Notably, these models can capture subtle features of folding pathways, such as the stabilization of misfolded states with specific structural characteristics. For chignolin, the machine-learned CG model successfully stabilized the same misfolded state with misaligned tyrosine residues as found in reference atomistic simulations [5]. This level of accuracy, combined with the computational efficiency that makes these models several orders of magnitude faster than all-atom simulations, positions machine-learned Cα models as powerful tools for exploring protein folding mechanisms.

Intrinsically Disordered Protein Dynamics

Coarse-grained models with implicit solvent have proven particularly valuable for studying intrinsically disordered proteins (IDPs), which challenge conventional experimental techniques due to their large-scale conformational fluctuations and transient structural elements [22]. The application of MARTINI 3 and SIRAH force fields to amyloid-β (Aβ42) monomers, with and without disulfide bonds, has revealed how structural constraints influence IDP behavior [22].

These studies demonstrated that the addition of a disulfide bond decreases the β-content of Aβ42 while increasing the tendency of the monomer to form fibril-like conformations, explaining various aggregation rates observed in experiments [22]. Analysis of monomer compactness, secondary structure content, and comparison between calculated and experimental chemical shifts confirmed that CG approaches provide reasonable methods to study IDPs, though they may lack some atomistic details due to the inherent simplifications [22].

Large-Scale Biomolecular Assemblies

Recent advances in CG-MD methodologies have enabled the simulation of massive biomolecular assemblies approaching cellular scales. The development of dynamic load balancing schemes in software like GENESIS CGDYN allows efficient simulation of non-uniform densities in large biological systems [21]. This capability has proven essential for studying biomolecular condensation through liquid-liquid phase separation (LLPS), where systems exhibit dramatic changes in density distribution over time.

Applications include simulations of IDP droplets, where researchers have observed fundamental processes such as droplet fusion and Ostwald ripening - a phenomenon where small droplets dissolve and their molecules redeposit into larger droplets [21]. These simulations achieve system sizes comparable to those observed in microscopy, providing unprecedented insight into mesoscopic biological processes at the residue-level description [21]. The ability to simulate such large systems while maintaining chemical specificity represents a significant milestone in bridging molecular and cellular scales in computational biology.

In the simulation of biomolecular systems, particularly for studying protein folding, coarse-grained (CG) molecular dynamics provides a crucial bridge between the computationally expensive detail of all-atom models and the biologically relevant time and length scales of functional processes [16] [6]. The core of developing any CG model lies in its parameterization—the derivation of effective forces that govern the interactions between CG sites. This process generally follows one of two fundamental philosophies: top-down or bottom-up parameterization [24] [25]. These approaches differ not only in their methodological foundations but also in their applications, strengths, and limitations. Top-down parameterization aims to reproduce experimental or macroscopic data, while bottom-up methods seek to derive parameters systematically from underlying, higher-resolution models [24]. For researchers focused on small protein folding, the choice between these philosophies significantly influences the model's transferability, computational efficiency, and physical accuracy. This article delineates these two parameterization strategies, their theoretical underpinnings, and provides practical protocols for their application in protein folding studies.

Theoretical Foundations

The Bottom-Up Approach

The bottom-up approach is grounded in statistical mechanics and aims to create a CG model whose properties are consistent with those of a more detailed, fine-grained (FG) reference system, typically an all-atom model [16] [24]. The central objective is to define an effective Hamiltonian for the CG model that faithfully represents the FG system by integrating out the unnecessary degrees of freedom.

The cornerstone of this approach is the potential of mean force (PMF), which serves as the physical basis for the coarse-grained force fields [16]. The PMF, denoted as ( W(\mathbf{R}) ), for a set of CG coordinates ( \mathbf{R} ) is defined by: [ W(\mathbf{R}) = -kB T \ln \left[ \frac{\int \exp\left(-V(\mathbf{R}, \mathbf{r}) / kB T\right) d\mathbf{r}}{\int d\mathbf{r}} \right] ] where ( V(\mathbf{R}, \mathbf{r}) ) is the potential energy of the atomistic system, ( \mathbf{r} ) represents the integrated-out atomistic degrees of freedom, ( k_B ) is Boltzmann's constant, and ( T ) is the absolute temperature [16]. The mean force acting on a CG site is then given by the negative gradient of the PMF, ( -\nabla W(\mathbf{R}) ).

In practice, the motion of CG sites is governed by a Langevin equation, which includes the mean force, a friction force, and a stochastic force: [ \frac{d}{dt}\mathbf{P}i(t) = -\nabla{\mathbf{R}i}W(\mathbf{R}) - \Gamma \dot{\mathbf{R}}i + \delta \mathbf{F}i^Q(t) ] The friction ( \Gamma ) and stochastic force ( \delta \mathbf{F}i^Q ) arise from the integrated degrees of freedom [16]. Several systematic methods exist to approximate the PMF and derive the CG force field:

• Variational Force Matching: This method trains a CG potential, ( U(\mathbf{x}; \boldsymbol{\theta}) ), by minimizing the difference between the CG forces and the mapped all-atom forces. The loss function is: [ L(\boldsymbol{\theta}) = \frac{1}{3nM} \sum{c=1}^{M} \| \boldsymbol{\Xi} \mathbf{F}(\mathbf{r}c) + \nabla U(\boldsymbol{\Xi}\mathbf{r}c; \boldsymbol{\theta}) \|^2 ] where ( \boldsymbol{\Xi} ) is the mapping operator, ( \mathbf{F}(\mathbf{r}c) ) are the all-atom forces, and the sum is over ( M ) configurations [6]. Machine learning potentials, particularly neural networks, are increasingly used to represent ( U ) [6].
• Relative Entropy Minimization: This approach minimizes the relative entropy, or Kullback-Leibler divergence, between the equilibrium probability distributions of the FG and CG systems [6].
• Inverse Boltzmann: One of the earliest methods, it determines pair potentials by inverting the radial distribution functions obtained from FG simulations [24].

The Top-Down Approach

In contrast, the top-down approach prioritizes the reproduction of experimental data or macroscopic properties over a strict connection to an atomistic force field [24] [25]. Parameters within a pre-defined analytical potential form are adjusted until the CG model's output matches target data, which can include thermodynamic properties, structural features, or known phenomenological behavior.

This philosophy is often employed in structure-based models, such as the Gō model, where the energy function is designed to stabilize the native structure of a protein [25]. Interactions present in the native state are assigned attractive potentials, while non-native contacts may be repulsive. Similarly, the Elastic Network Model (ENM) uses a simple harmonic potential to maintain the native structure: [ \frac{U}{kB T} = \sum{i \neq j} h(rc - r{ij}^0) \frac{\gamma}{2} (\Delta r)^2 ] where ( h ) is the Heaviside function, ( rc ) is a cutoff distance, ( r{ij}^0 ) is the native distance between atoms ( i ) and ( j ), and ( \gamma ) is the spring constant [7]. The parameters (e.g., ( \gamma ), ( r_c )) are not derived from atomistics but are tuned to reproduce experimental observables like B-factors or fluctuation profiles [7].

Knowledge-based statistical potentials also fall under the top-down umbrella. They are derived from databases of protein structures (e.g., the Protein Data Bank) by inverting frequency distributions of interactions: [ F{ij}(r) = -kT \ln \frac{P{ij}(r)}{P0} ] where ( P{ij}(r) ) is the observed probability of residue pair ( (i,j) ) being at distance ( r ), and ( P_0 ) is a reference probability [25].

Comparative Analysis: Strengths and Limitations

Table 1: Comparative analysis of bottom-up and top-down parameterization approaches.

Feature	Bottom-Up Approach	Top-Down Approach
Theoretical Basis	Statistical mechanics; Potential of Mean Force [16] [24]	Empirical fitting; Database statistics [25]
Primary Target	Consistency with underlying all-atom model forces and structure [6] [24]	Reproduction of experimental data (e.g., densities, yields) or specific structural features [24] [25]
Transferability	Higher in principle for states near the training reference; can extrapolate poorly to untrained states [6]	Often limited to the specific context or protein for which it was parameterized; can be system-specific [25]
Treatment of Dynamics	Can preserve correct thermodynamics by construction; dynamics require careful calibration of friction [16] [6]	Dynamics are often not physically meaningful; primarily used for structure sampling and thermodynamics [25]
Computational Cost	High (requires extensive all-atom simulation data for training) [6]	Low to moderate (relies on experimental data or simple potentials)
Best-Suited Applications	Studying processes where atomistic detail is crucial (e.g., ligand binding, ion effects) [24]; Multi-scale simulation	Rapid structure prediction; Studying large-scale conformational changes near the native state [7] [25]

The choice between bottom-up and top-down parameterization is dictated by the scientific question and available resources. The bottom-up approach excels when the research goal requires a direct, rigorous connection to the underlying atomistic physics. Its key strength is its foundation in statistical mechanics, which allows it, in principle, to preserve the correct thermodynamics of the system from which it was derived [6]. This makes it particularly suitable for studying processes where the influence of atomistic detail, such as specific ion effects or solvent-mediated interactions, is critical [24]. Furthermore, as machine learning coarse-grained potentials demonstrate, a single bottom-up potential can be trained on multiple proteins, showing promising transferability across different systems [6].

The primary limitations of the bottom-up approach are its computational cost—requiring extensive all-atom simulations for training—and the challenge of dynamics preservation. While thermodynamics can be captured by the PMF, the kinetic rates of processes are often accelerated, and friction terms must be carefully introduced to recover realistic dynamics [16].

The top-down approach, conversely, offers computational efficiency and ease of development. Models like the ENM are exceptionally powerful for studying collective motions and large-scale conformational changes near a protein's native state [7]. They are invaluable when the research focus is on deriving insights that are directly comparable to experimental observables, such as B-factors from crystallography.

The main drawback of top-down methods is their limited transferability. A model parameterized for a specific protein or under specific conditions may perform poorly outside that context [25]. Furthermore, because the energy landscape is not derived from a physical FG model, the energetics and dynamics may lack the granularity needed to study processes like folding pathways or ligand binding in detail. Knowledge-based potentials also face theoretical criticism regarding the physical meaning of their terms and the handling of many-body correlations [25].

Application Notes and Protocols

Protocol 1: Bottom-Up Force Matching for a Small Protein

This protocol outlines the development of a bottom-up CG model for a small protein using a force-matching approach, suitable for folding studies [6].

Workflow Diagram: Bottom-Up Force Matching

Step 1: Generate All-Atom Reference Data

• Procedure: Perform extensive, unbiased all-atom molecular dynamics simulations of the target protein(s) in explicit solvent. For folding studies, this should include simulations starting from both folded and unfolded states to ensure adequate sampling of the free energy landscape. The total simulation time may need to reach the microsecond-to-millisecond scale to observe multiple folding/unfolding events, as was done for a dataset of 12 proteins totaling ~9 ms [6].
• Deliverable: A trajectory file containing atomic coordinates and the corresponding forces for a large number of frames (e.g., millions of configurations).

Step 2: Define CG Mapping

• Procedure: Choose a mapping scheme that reduces the number of degrees of freedom. A common choice for proteins is to represent each amino acid by a single bead located at the Cα atom position [6] [25]. More complex mappings can include backbone and side-chain beads.
• Deliverable: A mathematical mapping operator ( \boldsymbol{\Xi} ) that transforms all-atom coordinates ( \mathbf{r} ) to CG coordinates ( \mathbf{x} = \boldsymbol{\Xi} \mathbf{r} ).

Step 3: Construct the CG Force Field

• Procedure: A neural network potential (NNP) is a powerful tool for representing the CG force field. The architecture should take the CG coordinates ( \mathbf{x} ) and particle types as input and output the total potential energy, ( U(\mathbf{x}; \boldsymbol{\theta}) ). The forces are then obtained by automatic differentiation: ( \mathbf{F}_{CG} = -\nabla U ).
• Prior Potentials: To improve learning and stability, incorporate simple analytical prior potentials for bonded terms (bonds, angles) and repulsive non-bonded terms to prevent chain rupture and atom clashes [6]. A chirality prior can also be enforced to prevent mirror-image structures [6].

Step 4: Train the Model via Force Matching

• Procedure: Minimize the force-matching loss function ( L(\boldsymbol{\theta}) ) (defined in Section 2.1) using a stochastic gradient descent algorithm. The training involves:
  • Feeding a batch of mapped CG coordinates ( \mathbf{x}c ) into the NNP.
  • Calculating the predicted CG forces ( -\nabla U(\mathbf{x}c; \boldsymbol{\theta}) ).
  • Comparing them to the reference mapped all-atom forces ( \boldsymbol{\Xi} \mathbf{F}(\mathbf{r}_c) ).
  • Updating the parameters ( \boldsymbol{\theta} ) to reduce the loss.
• Deliverable: A trained NNP that can predict energies and forces for any configuration of the CG system.

Step 5: Validation

• Procedure: Run CG-MD simulations using the trained potential. Validate the model by comparing its performance against data not used in training. Key metrics for folding studies include:
  • The population of the native state and other metastable states.
  • The free energy landscape as a function of relevant reaction coordinates (e.g., RMSD, fraction of native contacts).
  • Structural properties like radius of gyration.
• The model is successful if the thermodynamics match those of the all-atom reference system [6].

Protocol 2: Top-Down Structure-Based Model for Folding

This protocol details the creation of a top-down Gō-like model for studying the folding pathway of a small protein whose native structure is known [25].

Workflow Diagram: Top-Down Structure-Based Model

Step 1: Obtain the Native Structure

• Procedure: Acquire the high-resolution native structure of the target protein from the Protein Data Bank (PDB). This structure will serve as the single reference for defining the attractive interactions in the model.

Step 2: Identify Native Contacts

• Procedure: Analyze the PDB file to identify all pairs of CG sites (e.g., Cα atoms) that are considered to be in "native contact." A common criterion is a cutoff distance (e.g., 8-10 Å) between any two heavy atoms of different residues. This list of contacts defines which non-bonded interactions will be attractive.

Step 3: Define the Energy Function

• Procedure: Construct a Hamiltonian that favors the formation of native contacts and the native geometry. A typical form includes:
  • Bonded potentials: Harmonic potentials between consecutive Cα atoms, with equilibrium lengths set to their distances in the native structure.
  • Angular potentials: Harmonic potentials for angles defined by three consecutive Cα atoms, with equilibrium values from the native structure.
  • Non-bonded potentials: A Lennard-Jones-like potential for native contacts, with a minimum at the native distance. Non-native contacts are typically modeled with a purely repulsive potential.
• Deliverable: A complete analytical energy function parameterized by the native structure.

Step 4: Simulation and Refinement

• Procedure: Perform molecular dynamics or Monte Carlo simulations using the defined energy function. The initial simulation can start from a random coil or an extended structure.
• Parameter Tuning: If the simulation fails to fold to the native state or folds too easily, adjust the relative strengths of the energy terms (e.g., the well depth of the native contact potential). The model can be refined by comparing computed properties (e.g., stability, B-factors) with experimental data.

The Scientist's Toolkit: Key Research Reagents and Materials

Table 2: Essential software tools and resources for coarse-grained modeling of protein folding.

Tool/Resource	Type	Primary Function	Relevance to Approach
All-Atom MD Software (GROMACS, NAMD, AMBER)	Software Suite	Running high-resolution simulations to generate reference data for bottom-up training.	Bottom-Up
CafeMol	Simulator	A platform for simulating CG models, notably using Gō-like and other structure-based potentials.	Top-Down
AWSEM	CG Force Field	A physics-based CG force field that uses a knowledge-based potential; can be considered a hybrid approach.	Both
MARTINI	CG Force Field	A popular top-down force field parameterized to reproduce thermodynamic data; widely used for biomolecules in a membrane environment.	Top-Down
UNRES	CG Force Field	A physics-based CG model developed for protein folding and structure prediction.	Bottom-Up
Neural Network Libraries (PyTorch, TensorFlow)	Library	Building and training machine learning potentials for force matching and other bottom-up methods.	Bottom-Up
Protein Data Bank (PDB)	Database	Source of native protein structures for defining interactions in top-down, structure-based models.	Top-Down

Both bottom-up and top-down parameterization philosophies offer distinct and valuable pathways for constructing coarse-grained models of proteins. The bottom-up approach, with its rigorous foundation in statistical mechanics, is ideal for investigations where thermodynamic consistency with an all-atom reference is paramount, and where computational resources for generating training data are available. The emergence of machine-learned potentials is significantly advancing the power and accuracy of this approach [6]. Conversely, the top-down approach provides a computationally efficient and intuitive framework for exploring large-scale conformational dynamics and for generating structural hypotheses, especially when the native state is known but atomic detail is less critical.

For the researcher studying small protein folding, the choice is not necessarily binary. A promising future direction lies in hybrid approaches that leverage the strengths of both philosophies. A model could, for instance, use a bottom-up derived potential for the protein backbone while employing top-down informed potentials for specific long-range interactions. By understanding the core principles, strengths, and limitations of each method, scientists can make informed decisions to apply or combine these powerful parameterization strategies, thereby accelerating the discovery of protein folding mechanisms and their implications in health and disease.

Machine Learning Breakthroughs: Implementing Next-Generation CG Force Fields

Molecular dynamics (MD) simulations are a powerful tool for studying biomolecular processes, such as protein folding and ligand binding. However, the computational cost of all-atom (AA) simulations severely limits the time and length scales that can be accessed. Coarse-grained (CG) models, which reduce the number of degrees of freedom by grouping atoms into effective interaction sites, offer a solution by dramatically extending the scope of feasible simulations [5] [26]. Traditional CG models often rely on simplified, pairwise interaction potentials, which can fail to capture crucial many-body effects, leading to inaccuracies in representing the underlying free energy landscapes [27].

The integration of machine learning (ML), particularly deep learning, has instigated a paradigm shift in the development of CG force fields. ML approaches can learn complex, multibody potential functions directly from high-resolution reference data, such as all-atom simulations [28]. This article focuses on CGnets and other neural network architectures designed to model many-body interactions, providing application notes and detailed protocols for their use in studying small protein folding. These models achieve a critical balance, offering the computational efficiency of CG modeling while striving for the accuracy of finer-grained models [5] [27].

Theoretical Foundations and Key Architectures

The Force-Matching Principle

Bottom-up CG force fields aim to reproduce the equilibrium distribution of the reference all-atom system. A theoretically grounded method to achieve this is force-matching, also known as the multiscale coarse-graining (MS-CG) method [29] [27]. The core objective is to learn a CG force field whose predictions, for a given CG configuration, match the instantaneous forces from the atomistic system projected onto the CG coordinates (the local mean force). The loss function is the mean-squared error between the predicted and reference forces [29] [27]. When a neural network is trained to minimize this error, it effectively learns a potential of mean force (PMF) that approximates the many-body CG free energy surface [27].

Neural Network Architectures for CG Force Fields

Table 1: Key Neural Network Architectures for Coarse-Grained Force Fields.

Architecture	Core Description	Key Features	Demonstrated Applications
CGnet [27]	A deep neural network that uses hand-selected internal molecular features as input.	Enforces physical invariances; can be regularized with a prior energy to avoid unphysical states.	Alanine dipeptide; folding/unfolding of Chignolin [27].
CGSchNet [29]	Hybrid architecture combining CGnet with a SchNet graph neural network for feature learning.	Learns its own molecular features via continuous-filter convolutions; inherently more transferable.	Small biomolecular systems and the miniprotein Chignolin [29].
Equivariant GNNs (e.g., NequIP, Allegro) [30]	Graph neural networks that build E(3)-equivariance directly into the network layers.	Extremely data-efficient; can produce functional models with very limited data.	CG water [30].
QuinNet [31]	An end-to-end graph neural network that explicitly represents interactions up to five-body terms.	High expressiveness for complex many-body interactions without excessive computational cost.	Small molecules (MD17/22) and the miniprotein Chignolin [31].

The architectures in Table 1 address the many-body interaction challenge differently. CGnet provides a foundational framework where the network learns corrections to a physically motivated prior energy [27]. CGSchNet advances this by automating feature learning, which improves transferability across molecular systems [29]. A significant breakthrough comes from Equivariant Neural Networks, which explicitly encode the physical symmetry that the energy should be invariant to rotation and translation, and forces should be equivariant. This built-in physical prior drastically reduces the amount of training data required [30]. Finally, architectures like QuinNet explicitly design the network to capture very high-order (e.g., five-body) interactions, which are often necessary for modeling complex molecular systems with high accuracy [31].

Figure 1: A generalized workflow for developing and deploying a machine-learned coarse-grained force field, from all-atom reference data to production CG molecular dynamics simulations.

Application Notes: Performance and Validation

Recent studies demonstrate the powerful capabilities of ML-based CG models. A landmark 2025 study showed that a transferable CG model trained on a diverse set of proteins could successfully predict the metastable states of folded, unfolded, and intermediate structures for new protein sequences with low similarity (<40%) to the training set [5]. The model predicted folding/unfolding transitions, fluctuations of intrinsically disordered proteins, and relative folding free energies of protein mutants, achieving several orders of magnitude speed-up compared to all-atom MD [5].

For specific proteins like Chignolin, classical CG models with few-body interaction terms often fail to reproduce the folding/unfolding landscape. In contrast, CGnet successfully captures all free energy minima, highlighting the critical importance of the multibody terms it learns [27]. Furthermore, the data efficiency problem—where early NNPs required massive amounts of training data—is being overcome. Equivariant networks have been shown to produce functional models for CG water using datasets as small as a single frame of reference data, a significant improvement over non-equivariant models [30].

Table 2: Quantitative Performance of ML-CG Models on Benchmark Systems.

System	Architecture	Key Performance Metric	Result	Reference
Chignolin	CGnet	Reproduces folding/unfolding landscape	Captures all free energy minima, unlike classical few-body models	[27]
Chignolin	CGSchNet	Required training data	180 μs of reference AA simulation	[30]
8-peptides & fast-folding proteins	CGSchNet	Transferability (sequence similarity)	Successful prediction on sequences with <40% similarity to training set	[5]
CG Water	Equivariant GNN (NequIP)	Data efficiency	Functional model produced from a single frame of reference data	[30]
Chignolin	QuinNet	Accuracy on complex systems	Surpasses leading models on MD22 and Chignolin benchmarks	[31]

A persistent criticism of neural network potentials is their "black-box" nature. However, tools from explainable AI (XAI), such as Layer-wise Relevance Propagation (LRP), are now being applied to decompose the network's energy prediction into physically meaningful n-body contributions [32]. For example, applying GNN-LRP to a CG model of the protein NTL9 allowed researchers to pinpoint stabilizing and destabilizing interactions in various metastable states and interpret mutation effects, aligning the model's reasoning with established physical knowledge [32].

Figure 2: High-level architecture of a Graph Neural Network (GNN) based force field and its interpretation using Explainable AI (XAI). The GNN processes coarse-grained coordinates to predict energies and forces. Techniques like GNN-LRP can then decompose the prediction into contributions from specific n-body interactions, peering inside the black box [32].

Detailed Experimental Protocols

Protocol: Building a Transferable CG Model for Protein Folding

This protocol outlines the key steps for developing a transferable, machine-learned CG model for studying small protein folding, based on the methodology presented in [5].

• Curate a Diverse Training Set:

  • Select a set of small proteins and peptides with diverse folded structures (e.g., alpha-helices, beta-sheets, mixed motifs).
  • Generate extensive all-atom explicit-solvent MD simulations for these proteins. Ensure sampling covers folded, unfolded, and intermediate states where possible.
  • Include dimers of mono- and dipeptides to better model local interaction potentials.

• Define the Coarse-Grained Mapping:

  • A common and effective mapping for proteins is to represent each amino acid by a single bead located at the Cα atom [5] [27]. This drastically reduces the number of particles compared to an all-atom model.

• Prepare the CG Training Dataset:

  • Apply the CG mapping (e.g., Cα positions) to every frame of the all-atom simulations to generate the input CG configurations, X.
  • For each CG configuration, compute the instantaneous coarse-grained forces. This is the projection of the all-atom forces onto the CG coordinates (ΞF(F(r)) in Eq. 3 of [29]), which serves as the target Y for supervised learning.

• Select and Configure the Neural Network Architecture:

  • Choose an architecture such as CGSchNet [5] [29] for its transferability, or an Equivariant GNN [30] for superior data efficiency.
  • The input to the network is the set of CG coordinates for a given configuration.
  • The output of the network is the predicted potential energy of the CG system, U(x; θ).
  • The loss function L is the mean-squared error between the negative gradient of the predicted energy and the target forces: L = ‖ -∇U(x; θ) - ΞF(F(r)) ‖².

• Train and Regularize the Model:

  • Use a standard optimizer like Adam to minimize the loss function.
  • Incorporate a prior energy term (as in regularized CGnets [27]) to penalize unphysical configurations (e.g., atom clashes). This improves the stability of downstream simulations.
  • Split the data into training and validation sets. Use the validation loss for early stopping to prevent overfitting.

• Validate the Model on Unseen Proteins:

  • Test the trained model on small, fast-folding proteins not present in the training or validation sets (e.g., CLN025, TRP-cage, BBA, Villin headpiece) [5].
  • Run parallel-tempering or long Langevin dynamics simulations with the learned force field to compute the free energy landscape as a function of reaction coordinates like the fraction of native contacts (Q) and Cα root-mean-square deviation (RMSD).
  • Compare the predicted landscapes and metastable states with reference all-atom simulations or experimental data.

Protocol: Performing Free Energy Calculation with a Trained CG Model

This protocol describes how to use a trained ML-CG force field to study the folding thermodynamics of a small protein.

• System Setup:

  • Initialize the protein structure in an extended or unfolded state, defined by its CG coordinates.
  • Set up a Langevin dynamics integrator with an appropriate friction coefficient and a time step suitable for the CG model (typically larger than an all-atom time step).

• Enhanced Sampling:

  • To ensure converged sampling of folding/unfolding transitions, employ an enhanced sampling method. The Parallel Tempering (Replica Exchange) MD method is highly effective [5].
  • Run multiple replicas of the system at different temperatures (e.g., a range from 300 K to 500 K). Schedule exchanges between replicas periodically based on the Metropolis criterion.

• Production Simulation:

  • Run a sufficiently long simulation (or set of replica simulations) to observe multiple folding and unfolding events. The computational speed of the CG model makes this feasible on timescales equivalent to milliseconds or seconds of all-atom time.

• Trajectory Analysis:

  • Free Energy Surface (FES): Combine the trajectories from all replicas using the weighted histogram analysis method (WHAM) or a similar technique. Project the FES onto meaningful collective variables, such as:
    • Fraction of Native Contacts (Q): Measures the progress towards the native fold.
    • Cα Root-Mean-Square Deviation (RMSD): Measures the structural deviation from the native state.
  • Identify States: The minima on the FES correspond to metastable states (folded, unfolded, intermediates). Analyze representative structures from each minimum.
  • Compare Fluctuations: For the folded state, calculate the Cα root-mean-square fluctuations (RMSF) per residue and compare them to those from a reference all-atom simulation to assess the model's accuracy in capturing native dynamics [5].

The Scientist's Toolkit

Table 3: Essential Research Reagents and Computational Tools.

Item Name	Function / Description	Relevance to Protocol
All-Atom Force Field (e.g., CHARMM, AMBER)	Provides the high-resolution reference data through MD simulations.	Used to generate the target forces and configurations for training the CG model (Steps 1.1, 1.3).
CG Mapping Scheme (e.g., Cα-only)	Defines how atoms are grouped into coarse-grained beads.	Fundamental for defining the CG representation of the system (Step 1.2).
Neural Network Architecture (e.g., CGSchNet, Equivariant GNN)	The machine learning model that acts as the CG force field.	Core engine that learns the potential of mean force (Step 1.4).
Force-Matching Loss Function	The objective function that trains the NN to reproduce the reference forces.	Guides the optimization process during model training (Step 1.4).
Prior Energy Term	A simple physical potential that prevents the sampling of unphysical states.	Critical for stabilizing molecular dynamics simulations with the NN force field (Step 1.5) [27].
Parallel Tempering (Replica Exchange)	Enhanced sampling method to overcome free energy barriers.	Essential for obtaining converged folding/unfolding thermodynamics in production simulations (Step 2.2) [5].
Explainable AI (XAI) Tool (e.g., GNN-LRP)	Decomposes the NN's energy prediction into n-body contributions.	Used to interpret the model, build trust, and extract scientific insight (e.g., identifying key stabilizing interactions) [32].

The force-matching principle, also known as variational force-matching or the multiscale coarse-graining (MS-CG) method, provides a systematic, bottom-up approach for deriving accurate coarse-grained (CG) molecular models from high-resolution all-atom simulation data [6]. The core objective is to create CG potentials that preserve the thermodynamic properties of the all-atom system while accelerating molecular dynamics (MD) simulations by several orders of magnitude, making biologically relevant timescales accessible for small protein folding studies [6]. This data-driven approach has recently been enhanced through machine learning, particularly artificial neural network potentials (NNPs), which can learn many-body atomic interactions and construct a potential of mean force over the CG system without requiring pre-defined analytical potential forms [6].

For researchers investigating small protein folding, force-matching addresses a fundamental challenge: all-atom MD simulations, while providing atomic resolution, often cannot access the microsecond to millisecond timescales on which many proteins fold due to computational constraints [6] [33]. By creating CG representations that accelerate dynamics by more than three orders of magnitude while preserving system thermodynamics, force-matching enables the study of complete folding pathways and the identification of relevant structural states in the ensemble [6].

Theoretical Foundation of Force-Matching

Mathematical Formalism

The force-matching method is grounded in statistical mechanics and seeks thermodynamic consistency between the CG and all-atom models [6]. The approach requires a dataset ( {\Bbb{D}} ) of ( M ) coordinate-force pairs obtained from all-atom MD simulations: conformations ( {{{{{{{{\bf{r}}}}}}}}}{c}\in {{\Bbb{R}}}^{3{N}{c}} ) and corresponding forces ( {{{{{{{\bf{F}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}{c})\in {{\Bbb{R}}}^{3{N}{c}} ), where ( N_c ) is the number of atoms [6].

A linear mapping operator ( \Xi ) reduces the dimensionality of the atomistic system: [ {{{{{{{\bf{x}}}}}}}} = {{{{{{{\boldsymbol{\Xi }}}}}}}}{{{{{{{\bf{r}}}}}}}} \in {{\Bbb{R}}}^{3n} ] where ( 3n ) represents the remaining degrees of freedom in the CG representation [6]. Common mapping strategies include selecting only Cα atoms or backbone coordinates for protein systems.

The core of the method involves finding a potential ( U({{{{{{{{\bf{x}}}}}}}}}{c};{{{{{{{\boldsymbol{\theta }}}}}}}}) ) parameterized by ( \theta ) that minimizes the loss function: [ L({{{{{{{\bf{R}}}}}}}};{{{{{{{\boldsymbol{\theta }}}}}}}}) = \frac{1}{3nM}\mathop{\sum }\limits{c=1}^{M}\parallel {{{{{{{\boldsymbol{\Xi }}}}}}}}{{{{{{{\bf{F}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}{c}) + \nabla U({{{{{{{\boldsymbol{\Xi }}}}}}}}{{{{{{{{\bf{r}}}}}}}}}{c};{{{{{{{\boldsymbol{\theta }}}}}}}}){\parallel}^{2} ] This loss function ensures that the forces from the CG potential ( -\nabla U ) match the mapped all-atom forces ( {{{{{{{\boldsymbol{\Xi }}}}}}}}{{{{{{{\bf{F}}}}}}}} ) across all conformations in the training dataset [6].

Machine Learning Enhancement

Recent implementations have employed neural networks to represent ( U({{{{{{{{\bf{x}}}}}}}}}_{c};{{{{{{{\boldsymbol{\theta }}}}}}}}) ), creating neural network potentials (NNPs) that can learn complex many-body interactions [6]. These machine learning approaches can integrate multiple proteins into a single potential, demonstrating remarkable transferability across different protein systems and even to mutated sequences not included in the training data [6].

Table 1: Key Concepts in Force-Matching Theory

Concept	Mathematical Representation	Physical Interpretation
Mapping	( \mathbf{x} = \Xi \mathbf{r} )	Reduces system dimensionality by grouping atoms into CG beads
Force Matching	( \Xi \mathbf{F}(\mathbf{r}c) = -\nabla U(\Xi \mathbf{r}c; \theta) )	CG forces should match mapped all-atom forces
Loss Function	( L = \frac{1}{3nM} \sum{c=1}^M | \Xi \mathbf{F}(\mathbf{r}c) + \nabla U(\Xi \mathbf{r}_c; \theta) |^2 )	Quantifies discrepancy between CG and reference forces
Prior Potential	Bonded, repulsive, and chiral terms	Incorporates physical constraints to improve learning

Application to Small Protein Folding

Dataset Requirements and Generation

Successful application of force-matching to small protein folding requires extensive all-atom MD datasets that capture the relevant conformational space. A recent benchmark study created a unique dataset of unbiased all-atom MD simulations totaling approximately 9 milliseconds for twelve different fast-folding proteins with diverse secondary structure elements [6]. These proteins ranged from 10 to 80 amino acids and included various structural motifs:

• α-helical proteins: Polyalanine, WALP16 [33]
• β-sheet structures: WW Domain [6]
• β-turn structures: Chignolin, 1E0Q [6] [33]
• Mixed αβ structures: BBA, NTL9, Protein G, Trp-cage [6] [33]

For each protein, simulations began from random coil conformations and proceeded until the native structure was reached, ensuring coverage of both folded and unfolded states [6]. This comprehensive sampling is crucial for generating CG potentials that accurately represent the protein's energy landscape.

Table 2: Performance Comparison of Force-Matched CG Models for Protein Folding

Model Type	Speed Acceleration	Thermodynamic Preservation	Transferability	Typical Application
All-Atom MD	1x (baseline)	High (reference)	Limited to simulated system	Detailed mechanism studies
Force-Matched CG (NNP)	>1000x [6]	High (comparable energetics) [6]	Single potential for multiple proteins [6]	Folding pathways, state identification
United Atom Model	10-100x [34]	Moderate (mirrors key properties) [34]	Specific to chemical family	Bitumen, polymers [34]
Constrained MD	Order of magnitude [33]	Varies with clustering scheme	System-specific	Hierarchical folding studies [33]

Coarse-Grained Representation and Prior Potentials

In protein folding applications, common CG representations include:

• Cα-based mapping: Each amino acid is represented by a single bead at the Cα position [7]
• Backbone mapping: Multiple beads per residue to preserve backbone geometry
• Essential dynamics CG: Groups of atoms that move collectively form dynamic domains [7]

To enhance physical realism and improve learning efficiency, prior potentials are incorporated to constrain the conformational space [6]. These typically include:

• Bonded potentials: Maintain chain connectivity and prevent rupture
• Repulsive potentials: Prevent bead overlapping and unphysical clashes
• Chiral potentials: Enforce correct stereochemistry and prevent mirror images [6]

These priors reduce the complexity of the force field learning problem and can be viewed as imposing physical biases from domain knowledge, with the NNP performing delta-learning between the all-atom forces and the prior forces [6].

Experimental Protocol: Force-Matching for Protein Folding

The following diagram illustrates the comprehensive workflow for applying force-matching to develop coarse-grained models for protein folding studies:

Step-by-Step Protocol

Step 1: Reference All-Atom MD Simulations

• Objective: Generate comprehensive training data covering folded, unfolded, and intermediate states
• Procedure:
  • Select target proteins with diverse secondary structures (α-helical, β-sheet, mixed motifs)
  • Perform unbiased all-atom MD simulations starting from extended/unfolded conformations
  • Use explicit or implicit solvent models (e.g., GB/SA)
  • Simulate until native state is reached and multiple folding/unfolding events are observed
  • Collect coordinate-force pairs at regular intervals (e.g., every 1-10 ps)
• Validation: Ensure simulations recapitulate experimental structural features

Step 2: Coarse-Grained Mapping

• Objective: Define the reduced representation of the protein
• Procedure:
  • Select mapping resolution (Cα-only, backbone, or dynamic domains)
  • For Cα mapping: retain only Cα atoms as CG beads
  • Define bead identities based on amino acid type
  • Implement mapping operator Ξ to project all-atom coordinates to CG space
• Considerations: Balance between resolution reduction and preservation of essential structural features

Step 3: Prior Potential Design

• Objective: Incorporate physical constraints to guide learning
• Procedure:
  • Implement bonded potentials between adjacent beads in sequence: [ U{\text{bond}} = kb(r - r0)^2 ]
  • Add repulsive non-bonded potentials to prevent clashes: [ U{\text{rep}} = \epsilon\left(\frac{\sigma}{r}\right)^{12} ]
  • Enforce chirality through dihedral potentials: [ U{\text{chiral}} = k\phi[1 + \cos(n\phi - \delta)] ]
  • Parameterize prior terms based on statistical analysis of all-atom data

Step 4: Neural Network Potential Training

• Objective: Learn the CG potential that matches reference forces
• Procedure:
  • Design neural network architecture (typically multilayer perceptron or graph neural network)
  • Input: CG coordinates and bead type embeddings
  • Output: CG potential energy and forces
  • Train by minimizing force-matching loss function: [ L(\boldsymbol{\theta}) = \frac{1}{3nM} \sum{c=1}^M \| \Xi \mathbf{F}(\mathbf{r}c) + \nabla U(\Xi \mathbf{r}_c; \boldsymbol{\theta}) \|^2 ]
  • Use optimization algorithm (e.g., Adam, L-BFGS) with mini-batching
  • Regularize to prevent overfitting
• Validation: Monitor loss on validation set separate from training data

Step 5: Model Validation and Production Simulations

• Objective: Ensure CG model reproduces key structural and thermodynamic properties
• Procedure:
  • Perform CG-MD simulations starting from unfolded states
  • Verify recovery of experimental native structures
  • Compare free energy landscapes with all-atom reference
  • Validate against experimental data (NMR, FRET, etc.) not used in training
  • Assess transferability by simulating mutant sequences
• Success Metrics:
  • Correct folded state topology (RMSD < 2-3 Å)
  • Comparable state populations and thermodynamics
  • Acceleration factor > 1000x compared to all-atom MD

Table 3: Essential Computational Tools for Force-Matching Protein Folding Studies

Resource Category	Specific Tools/Solutions	Function/Purpose	Key Features
All-Atom MD Engines	AMBER, CHARMM, NAMD, GROMACS	Generate reference simulation data	Specialized hardware integration (GPUs), enhanced sampling methods
Coarse-Graining Software	VOTCA, BOCS, DeepCG, TensorMol	Perform force-matching optimization	Neural network potentials, loss function optimization, analysis tools
Neural Network Libraries	PyTorch, TensorFlow, JAX	Implement and train NNPs	Automatic differentiation, GPU acceleration, modular architecture
Analysis Packages	MDTraj, MDAnalysis, PyEMMA	Validate CG models and analyze folding	RMSD calculations, clustering, free energy estimation, state identification
Specialized Hardware	GPU Clusters, Anton2 Supercomputer	Accelerate reference and CG simulations	Massive parallelism, specialized processors for MD

Case Study: Multi-Protein Neural Network Potential

A recent landmark study demonstrated the power of force-matching by creating a single neural network potential capable of simulating twelve different proteins with diverse structural motifs [6]. This approach leveraged a massive dataset of approximately 9 milliseconds of all-atom simulation data across multiple proteins:

• Training Data: Unbiased all-atom MD of 12 fast-folding proteins (10-80 residues)
• CG Representation: Cα and side chain centers of mass
• Architecture: Neural network with residue-type embeddings
• Performance: >1000x acceleration while preserving thermodynamics
• Transferability: Single potential captured folding of all twelve proteins and predicted structural features of mutated sequences

This case study highlights the significant advantage of machine learning force-matching approaches: the creation of transferable potentials that capture universal aspects of protein energy landscapes rather than being limited to individual protein systems.

Analysis of Folding Pathways and Mechanisms

Force-matched CG models enable detailed analysis of protein folding pathways that would be computationally prohibitive with all-atom MD. Key analytical approaches include:

State Identification and Characterization

• Conformational Clustering: Use algorithms (k-means, hierarchical) to identify metastable states
• Principal Component Analysis: Project high-dimensional trajectories onto essential degrees of freedom
• Free Energy Estimation: Calculate state populations and relative stabilities

Folding Mechanism Elucidation

• Pathway Analysis: Identify folding intermediates and transition states
• Contact Formation Analysis: Track native and non-native contact formation along trajectories
• Hierarchical Folding: Investigate zipping-and-assembly mechanisms proposed for proteins like Trp-cage [33]

The following diagram illustrates how force-matched CG models reveal protein folding mechanisms through enhanced sampling:

The force-matching principle represents a powerful paradigm for developing coarse-grained models of protein folding that combine computational efficiency with physical accuracy. By directly learning from all-atom simulation data, these approaches preserve the essential thermodynamics of the system while accelerating dynamics by several orders of magnitude. The integration of machine learning, particularly neural network potentials, has further enhanced the transferability and accuracy of these models, enabling a single potential to describe multiple proteins and even predict the behavior of mutated sequences.

For researchers studying small protein folding, force-matched CG models provide unprecedented access to folding timescales and mechanisms, facilitating the identification of metastable states, folding pathways, and kinetic intermediates. As all-atom simulation datasets continue to grow and machine learning architectures become more sophisticated, force-matching promises to deliver increasingly accurate and transferable models that will further our understanding of protein folding and function.

The application of coarse-grained (CG) models in molecular dynamics (MD) has long been pursued to overcome the extreme computational cost of all-atom simulations, which remains a significant bottleneck in studying protein folding and dynamics [5]. A central challenge in the field has been the development of a transferable potential—a model parameterized on a set of training proteins that can accurately simulate the dynamics and thermodynamics of new, unseen protein sequences. Traditional CG approaches often required system-specific parameterization, limiting their predictive power for novel sequences [7]. The advent of machine learning (ML), combined with bottom-up coarse-graining principles, has recently transformed this landscape, making it feasible to develop models that learn the underlying physical principles of protein interactions rather than merely memorizing specific sequences or structures [5] [6]. This Application Note details the protocols and quantitative benchmarks for these next-generation, sequence-independent CG models, providing a framework for their application in small protein folding studies.

Key Advances in Machine-Learned Transferable Potentials

Recent breakthroughs have demonstrated that neural network potentials (NNPs) trained on diverse all-atom simulation data can achieve unprecedented transferability. The core innovation involves using deep learning to model the many-body interactions essential for realistic protein thermodynamics, which are difficult to capture with classical pair-wise potentials [5].

Table 1: Key Performance Metrics of Machine-Learned Coarse-Grained Models

Model / Study	Training Data	Transferability Demonstration	Reported Accuracy vs. All-Atom MD	Computational Speedup
CGSchNet [5]	Diverse all-atom protein simulations	Successful simulation of proteins with <40% sequence similarity to training set	Predicts metastable folded/unfolded states; matches fluctuations of disordered proteins; relative folding free energies consistent with experiment.	Several orders of magnitude
Multi-Protein NNP [6]	~9 ms of all-atom MD for 12 proteins	Single potential integrated all 12 proteins; captured structural features of mutated proteins	Preserved thermodynamics of systems; identified relevant structural states with comparable energetics.	>3 orders of magnitude
aSAM (Generative Model) [35]	MD datasets (ATLAS, mdCATH)	Generated ensembles for unseen proteins and temperatures	Replicated local flexibility (Cα RMSF PCC: ~0.89); captured temperature-dependent ensemble properties.	High computational efficiency vs. MD

The table above summarizes three seminal approaches. These models share a common methodology: they are trained using a bottom-up approach, often based on the variational force-matching method [6]. This technique trains the CG model to match the forces derived from all-atom reference simulations, ensuring thermodynamic consistency between the CG and all-atom ensembles.

Protocol: Developing and Applying a Transferable CG Model

This protocol outlines the workflow for developing a machine-learned, transferable CG potential and applying it to study the folding of a novel small protein.

Workflow Diagram

Step-by-Step Methodology

Stage 1: Data Curation and Coarse-Graining

• Generate Training Data: Assemble a large and diverse dataset of unbiased all-atom MD simulations. A representative dataset might include 12 fast-folding proteins with varied secondary structures (α-helices, β-sheets, mixed αβ), amounting to ~9 milliseconds of aggregate simulation time [6]. This diversity is critical for transferability.
• Define CG Mapping: Reduce the atomic resolution of the system. A common and effective mapping is to represent each amino acid by a single bead located at the Cα atom [5] [7]. This captures the essential dynamics of the protein backbone.

Stage 2: Neural Network Potential Training

• Architecture Selection: Employ a deep neural network architecture, such as CGSchNet [5] or other graph-based networks [6], which are well-suited for modeling molecular interactions.
• Loss Function and Training: Train the NNP using the variational force-matching method. The loss function minimized is: ( L(\boldsymbol{\theta}) = \frac{1}{3nM} \sum{c=1}^{M} \| \boldsymbol{\Xi} \mathbf{F}(\mathbf{r}c) + \nabla U(\boldsymbol{\Xi}\mathbf{r}c; \boldsymbol{\theta}) \|^2 ) where ( \boldsymbol{\Xi} \mathbf{F}(\mathbf{r}c) ) are the coarse-grained forces from all-atom data, ( \nabla U ) are the forces predicted by the NNP with parameters ( \boldsymbol{\theta} ), ( n ) is the number of CG sites, and ( M ) is the number of training configurations [6].
• Incorporate Physical Priors: Augment the NNP with simple analytical potential terms for bonds, angles, and repulsive interactions to prevent chain rupture and atom clashes, and to enforce correct chirality. This "delta-learning" approach simplifies the learning problem for the network [6].

Stage 3: Validation and Production Simulation

• Validate on Held-Out Proteins: Test the trained model on proteins not included in the training set, preferably those with low (<40%) sequence similarity [5]. Key validation metrics include:
  • Free energy surfaces as a function of reaction coordinates like fraction of native contacts (Q) and Cα root-mean-square deviation (RMSD).
  • Root-mean-square fluctuations (RMSF) of the folded state.
  • Relative folding free energies of mutants, if experimental data is available [5].
• Simulate Novel Proteins: Apply the validated model to simulate a novel protein sequence. Use enhanced sampling techniques like parallel tempering to efficiently sample the folding landscape [5].

The Scientist's Toolkit: Research Reagent Solutions

Table 2: Essential Computational Tools and Datasets for Transferable CG Modeling

Resource Type	Name / Example	Function & Application
CG Software/Code	CGSchNet [5], Multi-Protein NNP [6]	Provides the core architecture for the neural network potential; performs CG molecular dynamics simulations.
Reference Datasets	Multi-millisecond AA MD dataset [6], ATLAS[cite:8], mdCATH [35]	Serves as the ground-truth data for training bottom-up CG models via force-matching.
All-Atom MD Engines	OpenMM, GROMACS, AMBER, NAMD	Generates the reference all-atom simulation data required for training the CG model.
Analysis Tools	MDTraj, MDAnalysis, PyEMMA	Used to compute validation metrics like RMSD, RMSF, free energy surfaces, and contact maps from simulation trajectories.

Anticipated Results and Output Analysis

When applied to a novel small protein, a successful transferable CG model should yield the following results, benchmarked against all-atom MD or experimental data where possible:

Table 3: Quantitative Output Analysis and Benchmarking

Output Metric	Description	Interpretation & Benchmark
Free Energy Landscape	Surface calculated as a function of Q (fraction of native contacts) and Cα RMSD [5].	A well-folded protein will show a deep global minimum at high Q and low RMSD. Compare the stability of folded vs. unfolded basins to reference data.
RMSF (Cα)	Per-residue root-mean-square fluctuation within the folded basin [5].	Should correlate with flexibility from all-atom MD (e.g., Pearson correlation >0.88 [35]) and experimental B-factors.
Folding Mechanism	Pathways and intermediate states observed during simulations from extended/unfolded starting conditions [5].	The model should successfully fold the novel protein to a native-like structure and may identify key metastable intermediates.
Folding Free Energy (ΔG)	Relative stability calculated from the population of folded/unfolded states or via alchemical methods [5].	For mutants, predicted relative ΔG should trend with experimental stability data.

The workflow should successfully fold the novel protein, as illustrated in the following logic of the simulation process:

The development of transferable, sequence-independent CG models represents a paradigm shift in computational biophysics. By leveraging machine learning and bottom-up coarse-graining principles, these models now offer a feasible path to universal, quantitatively accurate, and computationally efficient protein simulations [5]. The protocols detailed herein provide a roadmap for researchers to implement these powerful tools, enabling the exploration of protein folding landscapes and the prediction of thermodynamic properties for novel sequences at a fraction of the computational cost of all-atom methods. Future advancements will likely focus on improving the modeling of complex motifs, integrating environmental conditions like pH and temperature more explicitly, and streamlining the backmapping process to recover all-atom detail [36] [35].

Within the context of developing accurate coarse-grained (CG) molecular dynamics (MD) parameters for small protein folding studies, the refinement of bonded interactions presents a significant challenge. All-atom (AA) MD simulations provide high-resolution insights into protein dynamics but are often computationally prohibitive for studying folding events, which can occur on timescales from microseconds to milliseconds [1]. Coarse-grained MD (CGMD) simulations, which simplify molecular structures into representative beads, offer a solution by dramatically extending the accessible timescales and length scales of simulations [37] [26]. The Martini force field, particularly its newest iteration Martini3, has emerged as a widely adopted top-down framework for CGMD simulations of biomolecules [37] [38]. While Martini3 provides excellent generality across molecular classes, its standardized parameters often lack the specificity required for accurately modeling complex processes like protein folding, where precise conformational dynamics are critical [37] [39]. This application note details a Bayesian optimization (BO) framework for refining Martini3 bonded parameters, enabling enhanced accuracy in protein folding studies while maintaining computational efficiency.

Bayesian Optimization Fundamentals

Bayesian optimization is a surrogate model-based algorithm designed to find the global optimum of black-box functions that are expensive to evaluate [40]. In the context of force field parametrization, each simulation to evaluate a set of parameters is computationally costly, making BO an ideal strategy. The algorithm operates by building a probabilistic surrogate model, typically a Gaussian process (GP), of the objective function based on observed data points [40]. This surrogate model is then used to construct an acquisition function that balances exploration (sampling regions with high uncertainty) and exploitation (sampling regions with promising predicted values) to suggest the next parameter set to evaluate.

The GP surrogate model provides a posterior distribution for function values at any test point ( \boldsymbol{x}* ), characterized by a mean function ( \mun(\boldsymbol{x}*) ) and variance ( \sigma^2n(\boldsymbol{x}*) ) [40]: [ \begin{aligned} \mun(\boldsymbol{x}*) &= K(\boldsymbol{x}, \boldsymbol{X}_n)[K(\boldsymbol{X}_n, \boldsymbol{X}_n) + \sigma^2 I]^{-1}(\boldsymbol{y} - m(\boldsymbol{X}_n)) + m(\boldsymbol{x}_) \ \sigma^2n(\boldsymbol{x}) &= K(\boldsymbol{x}_, \boldsymbol{x}*) - K(\boldsymbol{x}, \boldsymbol{X}_n)[K(\boldsymbol{X}_n, \boldsymbol{X}_n) + \sigma^2 I]^{-1}K(\boldsymbol{X}_n, \boldsymbol{x}_) \end{aligned} ]

The acquisition function, often Expected Improvement (EI), guides the selection of new points to evaluate [40]: [ \alpha{EI}(\boldsymbol{x}) = (\mu_n(\boldsymbol{x}_) - y^{best})\Phi(z) + \sigman(\boldsymbol{x})\phi(z) ] where ( z = \frac{\mu_n(\boldsymbol{x}_) - y^{best}}{\sigman(\boldsymbol{x}*)} ), and ( \Phi(\cdot) ), ( \phi(\cdot) ) are the standard normal CDF and PDF, respectively.

Workflow for BO-Driven Refinement of Martini3 Topologies

The following workflow illustrates the complete Bayesian optimization procedure for refining Martini3 bonded parameters, highlighting the integration of target property calculation, optimization logic, and iterative refinement crucial for achieving accurate parameters for protein folding studies.

Parameter Selection and Objective Definition

For protein folding applications, the bonded parameters defining molecular topology must be carefully selected. The core parameter set typically includes [37]: [ \boldsymbol{\theta} = {b0, kb, \Phi, k\Phi} ] where ( b0 ) represents bond lengths, ( kb ) bond constants, ( \Phi ) angle magnitudes, and ( k\Phi ) angle constants. For proteins containing aromatic residues, an additional parameter ( c ) may be included to account for constant aromatic bonds preserving ring topology [37].

The objective function must be designed to capture essential aspects of protein folding behavior. For small protein folding studies, key target properties include:

• Radius of gyration (( R_g )): Reports on protein compactness and dimensions
• Secondary structure formation: Monitoring α-helix and β-sheet content
• Native contact formation: Tracking development of native-like interactions
• Root mean square deviation (RMSD): Measuring deviation from native structure

These properties can be combined into a multi-objective function: [ \mathcal{L}(\boldsymbol{\theta}) = \sum{i=1}^{N} wi \left( \frac{Pi^{CG}(\boldsymbol{\theta}) - Pi^{ref}}{Pi^{ref}} \right)^2 ] where ( Pi^{CG}(\boldsymbol{\theta}) ) are properties from CG simulations with parameters ( \boldsymbol{\theta} ), ( Pi^{ref} ) are reference values from AA simulations or experimental data, and ( wi ) are weights balancing the relative importance of each property [37] [41].

Experimental Protocol: BO-Driven Refinement for a Small Protein System

Initial Setup and Reference Data Generation

Materials and Software Requirements:

Table 1: Essential Research Reagents and Computational Tools

Category	Item	Specification/Purpose
Simulation Software	GROMACS [38]	MD simulation engine (version 2019.x or later)
Coarse-Grained Force Field	Martini 3 [38]	Base CG force field for initial parameters
Bayesian Optimization	Custom Python implementation	BO algorithm implementation with Gaussian processes
Reference Data Generation	All-atom simulation system	AA force field (e.g., OPLS-AA, AMBER) with explicit solvent [38]
Analysis Tools	MDAnalysis, GROMACS tools	Trajectory analysis and property calculation
System Preparation	CGbuilder tool [38]	Mapping atomistic structures to CG representation

Step 1: Generate All-Atom Reference Data

• System Preparation: Obtain or generate the initial structure of the target small protein (e.g., 1-100 residues). For the villin headpiece (35 residues) or WW domain (~40 residues), PDB structures are publicly available [1].
• Solvation and Equilibration: Solvate the protein in an appropriate buffer (e.g., 150 mM NaCl) using ~
• Production Simulation: Run extensive AA MD simulations (≥100 ns) under relevant thermodynamic conditions to obtain converged reference data. For folding studies, include both folded and partially unfolded states if possible [1].
• Property Calculation: From the AA trajectories, calculate reference values for key properties: ( R_g ), native contact fraction, secondary structure content, and intramolecular distance distributions.

Step 2: Establish Coarse-Grained Mapping and Initial Parameters

• Mapping Scheme: Define the CG mapping for the protein. For Martini3, follow established mapping rules [38]:
  • Represent 4-5 heavy atoms with regular (R) beads for linear fragments
  • Use tiny (T) beads for conjugated, planar structures (e.g., aromatic rings)
  • Maintain chemical symmetry and group integrity
• Initial Topology: Construct the initial Martini3 topology using available building blocks or the Martini3 small molecule parametrization protocol [38].
• Parameter Bounds: Define reasonable search bounds for each bonded parameter based on Martini3 defaults and physical considerations (typically ±20-50% of initial values).

Bayesian Optimization Implementation

Step 3: Configure Bayesian Optimization

• Initial Design: Sample an initial space-filling design (e.g., Latin Hypercube) with 10-20 points across the parameter space.
• Surrogate Model: Implement a Gaussian process with Matérn 5/2 kernel to model the objective function [40].
• Acquisition Function: Select Expected Improvement (EI) as the acquisition function to balance exploration and exploitation.

Step 4: Iterative Optimization Loop

• Parallel Evaluation: For each parameter set in the current batch, run CGMD simulations:
  • System: Single protein in solution (matching AA reference conditions)
  • Duration: Sufficient to converge target properties (typically 100-500 ns)
  • Conditions: NPT ensemble, physiological temperature and pressure
• Property Extraction: Calculate target properties from CG trajectories.
• Objective Evaluation: Compute objective function value comparing CG and AA properties.
• Model Update: Update the GP surrogate model with new evaluation results.
• Candidate Selection: Maximize acquisition function to select next parameter set(s) for evaluation.
• Convergence Check: Terminate when objective function improvement falls below threshold or evaluation budget is exhausted.

Table 2: Example BO Parameters for Protein Folding Application

Parameter	Setting	Rationale
Initial Points	20	Sufficient to build initial surrogate model
Iterations	100-200	Balance between computational cost and convergence
Batch Size	5-10	Parallel evaluation to accelerate optimization
Convergence Tolerance	1% improvement over 20 iterations	Ensure stable convergence
Acquisition Function	Expected Improvement (EI)	Balanced exploration-exploitation

Validation and Production Simulations

Step 5: Validate Optimized Parameters

• Transferability Test: Validate optimized parameters on larger systems (e.g., higher degrees of polymerization) or related protein sequences [37].
• Folding Kinetics: If available, compare folding rates and pathways with experimental data or more detailed simulations.
• Thermodynamic Validation: Verify that the optimized model reproduces correct folding thermodynamics (e.g., melting behavior).

Step 6: Production Simulations

Execute large-scale or multiple folding simulations with optimized parameters to study folding mechanisms, kinetics, and pathways.

Key Considerations for Protein Folding Applications

Challenges in Protein Folding Simulations

Protein folding simulations present unique challenges that must be addressed in the parameter optimization process:

• Timescales: Folding events can range from microseconds to seconds, requiring enhanced sampling or very long simulations [1]
• Rare Events: The transition between folded and unfolded states may involve high free energy barriers
• Structural Heterogeneity: The unfolded ensemble contains diverse conformations that must be properly represented
• Force Field Accuracy: Energy landscapes must correctly stabilize the native state while properly representing unfolded and intermediate states [1]

Relationship Between Parameters and Folding Behavior

The following diagram illustrates how optimized bonded parameters influence molecular properties and ultimately affect the simulation of protein folding processes, highlighting the critical connections that enable accurate folding studies.

Multi-Objective Optimization for Folding

For protein folding applications, a single objective function may be insufficient. Consider a multi-objective BO approach that simultaneously optimizes [42]:

• Structural Accuracy: RMSD to native structure, native contact fraction
• Thermodynamics: Free energy landscape, folding temperature
• Kinetics: Folding rates, transition path times

The resulting Pareto front reveals trade-offs between objectives and allows selection of parameters based on the specific research goals.

Bayesian optimization provides a powerful, efficient framework for refining Martini3 bonded parameters specifically tailored for protein folding studies. By iteratively adjusting parameters to match reference data from all-atom simulations or experiments, the optimized CG models achieve accuracy comparable to AA simulations while maintaining the computational efficiency of Martini3. The protocol outlined here enables researchers to develop specialized force fields for studying folding mechanisms of small proteins, bridging critical gaps in our multiscale modeling capabilities. As BO methodologies continue to advance, particularly in handling high-dimensional parameter spaces [41], this approach promises to further enhance the accuracy and applicability of coarse-grained simulations for complex biomolecular processes.

The development of machine-learned coarse-grained (CG) models represents a significant advancement in the simulation of protein folding and dynamics. These models aim to bridge the gap between the high accuracy but extreme computational cost of all-atom molecular dynamics (MD) and the efficiency but often limited transferability of traditional CG approaches. By leveraging deep learning and large, diverse training sets, it is now feasible to construct transferable bottom-up CG force fields that capture the essential physics of protein thermodynamics and kinetics, enabling extrapolative molecular dynamics on novel protein sequences [5] [6]. This protocol details the practical steps for creating such models, from the initial design of the training dataset to the final deployment in production-scale simulations.

The overall process for developing and deploying a machine-learned coarse-grained model for protein folding studies can be broken down into several key stages, from data generation through to production simulation and analysis. The following diagram illustrates this integrated workflow.

Workflow for Machine-Learned Coarse-Grained Model Development. This diagram outlines the key stages in creating and deploying a transferable coarse-grained model for protein folding studies, from initial data generation to final production analysis.

Training Set Design and Data Generation

A foundational step in building a transferable CG model is the creation of a comprehensive and unbiased training dataset. The quality and diversity of this dataset directly determine the model's predictive capabilities and transferability.

Composition of a Diverse Protein Training Set

Table 1: Exemplary Protein Training Set Composition

Protein	PDB Code	Length (residues)	Secondary Structure	Simulation Time (per protein)
Chignolin	2RVD	10	β-hairpin	~0.1 - 1 ms
Trp-Cage	2JOF	20	α-helical	~0.1 - 1 ms
BBA	1FME	28	α/β mixed	~0.1 - 1 ms
Villin Headpiece	1YRF	35	α-helical	~0.1 - 1 ms
WW Domain	1E0L	37	β-sheet	~0.1 - 1 ms
Protein A (B-domain)	1BDD	60	α-helical bundle	~0.1 - 1 ms
Engrailed Homeodomain	1ENH	54	α-helical bundle	~0.1 - 1 ms
alpha3D	2A3D	73	α-helical bundle	~0.1 - 1 ms
NTL9	2HBA	56	α/β mixed	~0.1 - 1 ms
Protein G	1MI0	56	α/β mixed	~0.1 - 1 ms

The training set should encompass a wide range of secondary structure elements, including α-helices, β-sheets, and mixed α/β motifs, as exemplified in Table 1 [5] [6]. It is critical to include proteins with low sequence similarity (e.g., <40%) to any target proteins intended for future prediction, ensuring the model learns general physical principles rather than memorizing specific sequences [5]. The total dataset should amount to multiple milliseconds of aggregate simulation time to ensure sufficient sampling of folded, unfolded, and intermediate states [6].

All-Atom Reference Simulation Protocol

To generate the training data, perform unbiased all-atom MD simulations starting from extended or random coil conformations. This allows the model to learn the entire folding pathway.

• Force Field: Use a modern, well-validated all-atom force field (e.g., AMBER, CHARMM).
• Solvation: Employ explicit solvent models (e.g., TIP3P water) in a periodic box.
• Sampling Enhancement: Utilize generalized-ensemble algorithms to enhance conformational sampling where necessary. Parallel Tempering (Replica Exchange) is highly effective, where the exchange probability between replicas i and j is given by: ( P(Ci \rightarrow Cj) = \min(1, \exp(-\betaiE(Cj) - \betajE(Ci) + \betaiE(Ci) + \betajE(Cj))) ) [43].
• Output: Extract and save trajectories at a sufficient frequency to capture relevant dynamics. The final dataset consists of a large number of coordinate-force pairs (( \mathbf{r}c, \mathbf{F}(\mathbf{r}c) )) for the all-atom system [6].

Coarse-Grained Model Development

With a robust training dataset in place, the next phase involves selecting a model architecture and a learning objective to derive the CG potential.

CG Mapping and Prior Energy Terms

A common and effective CG mapping, ( \mathbf{x} = \boldsymbol{\Xi}\mathbf{r} ), is to retain only the Cα atoms of the protein backbone, thereby dramatically reducing the number of degrees of freedom [6]. To stabilize the chain and prevent unphysical conformations, it is essential to incorporate analytical prior energy terms. These terms reduce the complexity of the force-matching problem and embed basic physical constraints.

• Bonded Potential: A harmonic potential between consecutive Cα atoms to maintain chain connectivity.
• Angle Potential: A harmonic potential on virtual Cα-Cα-Cα angles to guide local chain geometry.
• Repulsive Potential: A soft-sphere repulsive term between non-bonded beads to prevent atomic clashes.
• Chirality Potential: A dihedral prior term to enforce correct backbone chirality and prevent mirror-image structures [6].

The neural network potential then effectively performs "delta-learning," correcting and augmenting these simple priors with learned multi-body interactions.

Neural Network Potential and Training

The core of the model is a neural network that learns a potential energy function, ( U(\mathbf{x}; \boldsymbol{\theta}) ), for the CG system.

• Architecture: Graph Neural Networks (GNNs) or SchNet are particularly well-suited as they are naturally equivariant and can handle molecular structures [5] [6]. Recent work also shows promise with standard transformer blocks, challenging the need for complex, domain-specific architectures [44].
• Learning Objective: The network is trained using the variational force-matching method. The loss function ( L(\boldsymbol{\theta}) ) is minimized: ( L(\boldsymbol{\theta}) = \frac{1}{3nM} \sum{c=1}^{M} \| \boldsymbol{\Xi} \mathbf{F}(\mathbf{r}c) + \nabla U(\boldsymbol{\Xi}\mathbf{r}_c; \boldsymbol{\theta}) \|^2 ) [6]. This ensures that the forces from the CG potential (( -\nabla U )) match the mapped all-atom forces as closely as possible.
• Training: The model is trained on the dataset of mapped coordinates and forces. A single potential can be trained on multiple proteins simultaneously to enhance its transferability [6].

Model Validation and Benchmarking

Before deploying a model in production, rigorous validation against both all-atom simulations and experimental data is mandatory.

Quantitative Validation Metrics

Table 2: Key Metrics for CG Model Validation

Validation Category	Specific Metric	Target Performance	Validation Method
Structural Accuracy	Cα RMSD to Native State	< 2.0 Å (small proteins)	PT/MD Simulation
	Fraction of Native Contacts (Q)	> 0.8 (folded state)	PT/MD Simulation
Thermodynamic Accuracy	Folding Free Energy (ΔG)	Match all-atom/expt. within ~kBT	PT/MD Simulation (FES)
	Free Energy Surface (FES) Topography	Reproduce metastable states	FES Projection (e.g., Rg vs RMSD)
Kinetic Accuracy	Folding/Unfolding Rates	Order-of-magnitude agreement	Long Langevin MD
Transferability	Stability of Unseen Proteins	Correct folded state as global min.	Simulation from extended chain

Validation should assess the model's ability to reproduce the thermodynamics and kinetics of proteins not included in the training set.

• Free Energy Surfaces: Use Parallel Tempering or long Langevin dynamics simulations to compute converged Free Energy Surfaces (FES) as a function of collective variables like the root-mean-square deviation (RMSD) from the native structure or the fraction of native contacts (Q) [5]. A successful model will recapitulate the key metastable basins (folded, unfolded, intermediates) found in all-atom references.
• Structural and Fluctuation Analysis: Compare the Cα root-mean-square fluctuations (RMSF) of the folded state from CG simulations to those from all-atom simulations to ensure the model captures the native state's flexibility correctly [5].
• Extrapolation Tests: Simulate the folding of larger proteins (e.g., 1ENH, 2A3D) from extended initial conditions. The model should spontaneously fold these proteins to structures close to their native state, demonstrating true transferability [5].

Production Simulations

Once validated, the model can be deployed for production simulations to investigate novel biological questions.

Simulation Setup and Execution

• Initial Conditions: Production runs can start from experimentally determined structures (to study dynamics) or from extended chains (to study folding).
• Integrator: Use a Langevin dynamics integrator with a friction coefficient and noise term that satisfies the fluctuation-dissipation theorem. This provides a thermostat and accounts for the implicit solvent.
• Timestep: The CG potential allows for significantly larger timesteps compared to all-atom MD, typically in the range of 10-100 femtoseconds [26].
• Sampling: For equilibrium properties, long constant-temperature simulations are sufficient. To rapidly explore complex free energy landscapes with high barriers, Parallel Tempering (Replica Exchange) remains the method of choice [5] [43].
• Performance: The resulting CG model is typically several orders of magnitude faster than equivalent all-atom MD, enabling the simulation of processes like folding on timescales that are otherwise inaccessible [5].

Analysis of Production Data

Analysis of production trajectories focuses on extracting thermodynamic and mechanistic insights.

• State Identification: Use clustering algorithms on the trajectories to identify dominant conformational states.
• Pathway Analysis: Tools like Principal Component Analysis (PCA) can be applied to folding trajectories to reduce dimensionality and identify the essential collective coordinates that capture the main structural changes along the folding pathway [26].
• Free Energy Calculation: The free energy as a function of any collective variable z can be computed from the simulation trajectory as ( F(z) = -k_B T \ln P(z) ), where ( P(z) ) is the probability distribution of z [43].

The Scientist's Toolkit

Table 3: Essential Research Reagents and Computational Tools

Tool/Reagent	Category	Primary Function	Exemplars / Notes
All-Atom MD Software	Simulation Engine	Generate training data	GROMACS, AMBER, NAMD, OpenMM
CG Model Software	Simulation Engine	Run production CG simulations	Custom Python/TensorFlow code, OpenMM for inference
Neural Network Library	Model Framework	Define and train NN potential	TensorFlow, PyTorch, JAX
UNRES Force Field	Traditional CG Model	Physics-based CG benchmark	~4000x faster than all-atom MD [26]
Martini Force Field	Traditional CG Model	CG for intermolecular interactions	Limited for intramolecular dynamics [5]
AWSEM/UNRES	Traditional CG Model	CG for folding & dynamics	Can struggle with alternative states [5]
Rosetta	Modeling Suite	De novo design & structure prediction	Used for sampling shape space [45]
Parallel Tempering	Sampling Algorithm	Enhance sampling of rugged landscapes	Critical for converging FES [5] [43]
Principal Component Analysis (PCA)	Analysis Tool	Identify essential dynamics from trajectories	Reveals key folding modes [26]

The study of fast-folding proteins provides critical insights into the fundamental principles of protein folding. Among these, chignolin, Trp-cage, and villin headpiece have emerged as model systems due to their small sizes and rapid folding kinetics. These proteins serve as essential benchmarks for developing and validating computational methods, particularly coarse-grained molecular dynamics simulations. The advent of machine-learned coarse-grained models represents a significant advancement, enabling accurate simulation of protein folding landscapes that were previously inaccessible with all-atom models due to computational constraints. This application note details successful folding case studies and experimental protocols, providing researchers with practical frameworks for implementing these methods in drug development and protein engineering applications.

Structural Characteristics and Folding Performance

Table 1: Structural Properties and Folding Performance of Model Proteins

Protein	Residues	Native Structure	Folded State Accuracy (RMSD)	Key Stabilizing Features
Chignolin	10	β-hairpin/turn	<1.0 Å [5]	β-hairpin formation, native contacts
Trp-cage	20	Two α-helices	<1.0-2.0 Å [46] [47]	Pro-Trp stacking, Asp-9/Arg-16 salt bridge, proline/glycine enrichment
Villin Headpiece	36	Three α-helices	~0.5 Å [5]	Hydrophobic core, specific tertiary interactions

Quantitative Folding Results from Computational Studies

Table 2: Quantitative Folding Results from Computational Studies

Simulation Method	Chignolin	Trp-cage	Villin Headpiece	Computational Speed
All-Atom MD (Anton)	~μs-ms timescale [48]	~μs-ms timescale [48]	~μs-ms timescale [48]	Reference baseline
Accelerated MD (aMD)	Complete folding captured [48]	Complete folding captured [48]	Complete folding captured [48]	Significant acceleration over cMD
Machine-Learned CG	Folds correctly, captures misfolded state [5]	Folds correctly [5]	Folds correctly [5]	3+ orders of magnitude faster than all-atom [6]
Discrete MD (DMD)	N/A	<2.0 Å RMSD [46] [47]	N/A	High efficiency via simplified potential

Experimental Protocols

Machine-Learned Coarse-Grained Force Field Protocol

Diagram 1: Workflow for Developing Machine-Learned Coarse-Grained Force Fields

Training Data Generation

• Source Data: Generate all-atom molecular dynamics simulations of the target proteins using explicit solvent models. For comprehensive training, include multiple secondary structure types and various protein sequences [6].
• Data Volume: Successful models have been trained on approximately 9 milliseconds of cumulative all-atom simulation data across twelve different proteins [6].
• Diversity Requirement: Ensure training data includes folded, unfolded, and intermediate states to adequately capture the complete free energy landscape.

Coarse-Grained Mapping

• Standard Mapping: Employ Cα atoms or four beads per residue (N, Cα, C, and side-chain centroid) to represent each amino acid [49] [6].
• Resolution Selection: Balance computational efficiency with structural accuracy. Higher resolution preserves more structural details but reduces computational gains.

Neural Network Force Field Training

• Architecture: Implement deep neural networks (e.g., CGnets, CGSchNet) that respect physical invariances (rotational and translational invariance) [5] [27].
• Loss Function: Utilize force-matching loss: (L(\boldsymbol{\theta}) = \frac{1}{3nM}\sum{c=1}^{M}\parallel \boldsymbol{\Xi}\mathbf{F}(\mathbf{r}c) + \nabla U(\boldsymbol{\Xi}\mathbf{r}_c;\boldsymbol{\theta})\parallel^2) where (\boldsymbol{\Xi}) is the mapping operator, (\mathbf{F}) are all-atom forces, and (U) is the CG potential [6] [27].
• Prior Knowledge: Incorporate physical prior terms for bonded interactions, repulsive forces, and chirality to prevent unphysical structures [6].

Model Validation

• Thermodynamic Consistency: Validate by comparing free energy surfaces and metastable state populations with all-atom references [6] [27].
• Transferability Testing: Assess performance on protein sequences not included in the training set, ensuring the model generalizes beyond training data [5].

Accelerated Molecular Dynamics (aMD) Protocol

System Setup

• Force Field Selection: Choose appropriate force fields (AMBER, CHARMM, etc.) that have been validated for the target proteins. Different force fields may bias toward specific secondary structures [48].
• Solvent Model: Employ explicit solvent models (e.g., TIP3P) for accurate solvation effects, as implicit models may over-stabilize certain conformations [48] [50].

aMD Parameters

• Dual-Boost Application: Apply both dihedral and total energy boosts using parameters tuned to the system size [48]:
  • Dihedral boost: (E{\text{dihed}} = V{\text{dihedavg}} + a1 \times N{\text{res}}), (\alpha{\text{dihed}} = a2 \times N{\text{res}}/5)
  • Total boost: (E{\text{total}} = V{\text{totalavg}} + b1 \times N{\text{atoms}}), (\alpha{\text{total}} = b2 \times N{\text{atoms}})
• Reweighting: Use cumulant expansion to the 2nd order to recover canonical distributions from aMD simulations [48].

Analysis Methods

• Free Energy Calculation: Construct free energy profiles along reaction coordinates such as RMSD and fraction of native contacts [48].
• State Identification: Identify and characterize unfolded, intermediate, and folded states through cluster analysis [48].

The Scientist's Toolkit

Table 3: Key Research Reagent Solutions for Protein Folding Studies

Resource	Type	Function/Application	Examples/Notes
Force Fields	Software Parameter Set	Defines molecular interactions	AMBER ff14SBonlysc, CHARMM22, AMBER ff99SBILDN [48]
Coarse-Grained Mapping	Computational Method	Reduces system complexity	Cα-only, 4-bead per residue (N, Cα, C, side-chain) [49] [6]
Enhanced Sampling Algorithms	Computational Method	Accelerates rare events	aMD, Replica Exchange MD (REMD) [48] [50]
Neural Network Architectures	Software Tool	Learns CG force fields	CGnets, CGSchNet [5] [27]
Validation Metrics	Analytical Method	Assesses model accuracy	RMSD, fraction of native contacts, free energy profiles [48] [5]

The successful folding of chignolin, Trp-cage, and villin headpiece using coarse-grained molecular dynamics represents a significant milestone in computational biophysics. Machine-learned coarse-grained force fields have demonstrated remarkable capability in capturing protein folding thermodynamics while achieving computational speedups of three orders of magnitude compared to all-atom simulations. These advances provide researchers with powerful tools for studying protein folding mechanisms, predicting mutant effects, and potentially designing novel protein structures. As these methods continue to evolve, they promise to bridge the gap between computational efficiency and physical accuracy, opening new possibilities for protein engineering and drug development applications.

Parameterization Pitfalls: Solving Common Challenges in CG-MD Simulations

The study of protein folding is a fundamental challenge in structural biology, with profound implications for understanding cellular function and developing therapeutic interventions. Coarse-grained (CG) molecular dynamics has emerged as a powerful strategy to access the biologically relevant timescales of folding events that remain beyond the reach of all-atom simulations. By reducing molecular complexity, CG models extend simulation capabilities to larger systems and longer timescales; however, they introduce significant sampling limitations due to the simplified energy landscapes that can trap simulations in non-physiological states. The rugged energy landscape of biomolecules features numerous local minima separated by high-energy barriers, making adequate sampling of conformational space a persistent challenge [51] [52]. This application note examines enhanced sampling techniques specifically adapted for CG landscapes within the context of small protein folding studies, providing researchers with practical protocols to overcome these limitations and achieve more accurate thermodynamic and kinetic characterization of folding processes.

The theoretical foundation of these methods rests upon the energy landscape theory of protein folding, which posits that biomolecular dynamics are governed by complex, hierarchical landscapes where conformational substates are separated by barriers that can exceed 100 kJ/mol [51]. Experimental studies using techniques such as fluorescence spectroscopy, nuclear magnetic resonance, and single-molecule experiments confirm the wide range of timescales involved in protein dynamics, from local atomic fluctuations to large-scale conformational changes [51]. Enhanced sampling methods provide the computational means to bridge the gap between these experimental observables and atomic-level behavior by facilitating barrier crossing in silico.

Enhanced Sampling Methodologies for CG Landscapes

Monte Carlo-Based Approaches

Monte Carlo methods provide a flexible framework for enhanced sampling in CG simulations, particularly valuable for their ability to handle discontinuous potentials and complex boundary conditions. The fundamental principle underlying canonical MC sampling is the generation of conformational ensembles consistent with Boltzmann statistics through the Metropolis acceptance criterion. This criterion dictates that states decreasing the system energy are always accepted, while those increasing energy are accepted with probability P = exp(-βΔU), where β = 1/kBT and ΔU is the energy difference between new and old configurations [51]. This controlled acceptance of energetically unfavorable moves enables the system to overcome barriers in the vast conformational space and escape local minima that would trap conventional simulations.

Several specialized MC variants have been developed to address specific challenges in CG protein folding studies:

• Configurational Bias MC (CB-MC): This approach enhances sampling efficiency for chain-like molecules by preferentially "growing" molecular segments toward favorable states. Unlike traditional random trial moves, CB-MC employs a Rosenbluth acceptance criterion where P = min(1, Wnew/Wold), with the Rosenbluth factor W calculated as the product of Boltzmann weights for trial segment insertions [51]. For certain applications such as chromatin folding, specialized implementations like "end-transfer CB-MC" can achieve dramatic efficiency improvements—quadratic versus exponential scaling with polymer length [51].

• Parallel Tempering (Replica Exchange): This method employs multiple simultaneous simulations at different temperatures, periodically exchanging configurations between replicas with a probability that maintains each temperature's equilibrium distribution: P = min{1, exp[(βnew - βold)(Unew - Uold)]}, where βi = 1/(kB Ti) and Ui is the internal energy [51]. The high-temperature replicas can cross energy barriers more easily, effectively feeding diverse configurations to the lower-temperature simulations. This approach has proven particularly effective in its molecular dynamics incarnation, termed replica exchange MD [51] [52].

• Wang-Landau Sampling: This algorithm performs multiple random walks in energy space, each sampling a different energy range, with the combined information enabling calculation of canonical averages and thermodynamic quantities at any temperature [51]. When compared with parallel tempering for protein conformational sampling, both methods demonstrated similar performance and were two orders of magnitude faster than canonical MC simulations at low temperature [51].

Table 1: Comparison of Monte Carlo Sampling Techniques for CG Protein Folding

Method	Key Mechanism	Acceptance Criterion	Advantages	System Suitability
Canonical MC	Metropolis criterion with random moves	P = min{1, exp(-βΔU)}	Simple implementation; general applicability	Small systems; discontinuous potentials
Configurational Bias MC	Biased growing of molecular segments	P = min(1, Wnew/Wold)	Efficient for chain molecules; reduces overlap	Polymers; chromatin; protein folding
Parallel Tempering	Temperature-based replica exchange	P = min{1, exp[(βnew-βold)(Unew-Uold)]}	Effective barrier crossing; natural parallelism	Rough energy landscapes; small proteins
Wang-Landau	Multiple random walks in energy space	Density of states-based	Direct thermodynamic calculation; single-processor efficiency	Thermodynamic property calculation

Modern implementations have successfully integrated MC modules into standard simulation packages like CHARMM, with automated optimization of step sizes and efficient combination with minimization or MD modules [51]. These optimized MC methods have demonstrated superior performance compared to standard Langevin dynamics in reaching folded states of small proteins, arguing for their continued development and application in biomolecular studies.

Harmonic Approximation Techniques

Harmonic approximation methods provide an efficient alternative for exploring protein dynamics by leveraging the intrinsic structural properties of biomolecules to reduce computational complexity. These approaches are particularly valuable for identifying collective motions relevant to folding processes and functional dynamics without the expense of full conformational sampling.

• Normal Mode Analysis (NMA) and Elastic Network Models: These techniques utilize spectral decomposition of a mass-weighted Hessian matrix at thermal equilibrium to identify natural vibration modes of a protein structure [51]. Modern extensions like elastic network models forgo computationally demanding diagonalization through simplified bead-spring representations that assume the native structure represents a minimum energy state [51] [7]. The potential energy in a standard elastic network is typically expressed as U/kBT = Σi≠j h(rc - rij0)(γ/2(Δr)2), where h(r) is the Heaviside function, Δr = rij - rij0, rij0 is the native distance, and γ is a spring constant [7]. These models have proven remarkably effective in capturing large-scale conformational changes relevant to folding using only Cα coordinates from native structures.

• Principal Component Analysis (PCA) and Essential Dynamics (ED): These methods focus on low-frequency, high-amplitude vibrational modes that often correspond to functionally relevant motions. Rather than assuming thermal equilibrium, PCA constructs a variance-covariance matrix of positional fluctuations from simulation trajectories, projecting configurations onto principal components that capture the most significant collective motions [51]. Applications have revealed how seemingly similar structures can exhibit different dynamic preferences—for example, PCA of TATA-box DNA sequences bound to TATA-binding protein demonstrated that high-efficiency variants favored complexation motions while low-efficiency variants tended toward dissociation deformations, despite remarkably similar structures [51].

• Dynamic Elastic Network Model (DENM): This extension of traditional ENM specifically addresses protein unfolding processes by incorporating force-dependent bond breaking and formation probabilities. The model describes noncovalent bond breaking under force loading with probability Poff = 1 - exp(-koffcosθΔtexp(γ(Δr)2/kBT)), while broken bonds may reform with probability Pon = 1 - exp(-koncosθΔtexp(-γ(Δr)2/kBT)) [7]. This approach successfully reproduces force-extension curves observed in atomic force microscopy experiments and all-atom MD simulations, providing a computationally efficient framework for studying mechanical unfolding.

Table 2: Harmonic Approximation Methods for CG Simulations

Method	Theoretical Foundation	Key Applications	Computational Cost	Limitations
Normal Mode Analysis (NMA)	Hessian matrix diagonalization	Collective motions; vibration modes	Moderate (diagonalization)	Harmonic assumption; single minimum
Elastic Network Model (ENM)	Simplified bead-spring potential	Large-scale motions; mechanical properties	Low (cutoff-based)	Native-centric; simplified interactions
Principal Component Analysis (PCA)	Dimensionality reduction of fluctuations	Essential dynamics; conformational ensembles	Low to moderate (covariance)	Sampling-dependent; trajectory quality
Dynamic ENM (DENM)	ENM with kinetic bond breaking	Unfolding pathways; force spectroscopy	Moderate (event monitoring)	Parameter sensitivity; Bell model assumptions

Machine Learning-Enhanced Sampling

The integration of machine learning with enhanced sampling represents a paradigm shift in coarse-grained modeling, addressing fundamental limitations in both sampling efficiency and potential accuracy. Recent approaches leverage neural network potentials (NNPs) to capture many-body interactions while maintaining computational efficiency sufficient for folding timescales.

A primary challenge in bottom-up coarse-graining is achieving thermodynamic consistency—ensuring the equilibrium distribution sampled by the CG model matches that of the all-atom reference system. The force matching approach addresses this by minimizing the loss function L(R;θ) = (1/3nM) Σc=1M ∥ΞF(rc) + ∇U(Ξrc;θ)∥², where Ξ is the mapping operator, F(rc) are reference forces, and U(Ξrc;θ) is the CG potential parameterized by θ [6]. This method allows a single potential to integrate multiple proteins while capturing experimental structural features, including those of mutated sequences not included in training [6].

A particularly innovative application of machine learning addresses the critical challenge of mapping entropy calculation, which quantifies the information loss associated with a specific coarse-graining scheme. Traditional computation of mapping entropy is computationally intensive, but deep graph networks (DGNs) can accelerate this calculation by a factor as large as 10⁵ while maintaining accuracy [53]. This dramatic speedup enables reconstruction of the complete mapping entropy landscape using Wang-Landau sampling, facilitating identification of optimal CG representations that preserve the maximal information about system behavior [53].

Enhanced sampling also improves the training process itself. By biasing along CG degrees of freedom during data generation, then recomputing forces with respect to the unbiased potential, researchers can simultaneously shorten the simulation time needed to produce equilibrated data and enrich sampling in transition regions while preserving the correct potential of mean force [54]. This strategy has demonstrated notable improvements for both model systems like the Müller-Brown potential and biological molecules such as capped alanine [54].

Experimental Protocols

Protocol 1: Replica Exchange for CG Protein Folding

Objective: Enhance sampling of folding landscapes for small proteins (10-80 amino acids) using temperature-based replica exchange.

Materials and Setup:

• CG Representation: Cα-based mapping with one bead per amino acid
• Prior Potentials: Include bonded terms (bonds, angles), repulsive interactions, and chirality constraints to prevent chain rupture and mirror image exploration [6]
• Temperature Distribution: 16-64 replicas exponentially spaced between 300K and 500K
• Simulation Package: GROMACS, NAMD, or CHARMM with REMD capabilities

Procedure:

• System Preparation:
  • Obtain native structure from Protein Data Bank
  • Generate CG coordinates using Cα atoms
  • Define prior potentials: bonds between sequential residues, angles, and repulsive non-bonded interactions
  • Apply chirality constraints using dihedral terms to prevent mirror image conformations [6]

• Replica Initialization:

  • For each replica, assign starting temperatures according to exponential distribution
  • Initialize structures from extended chain or native state (for folding/unfolding studies respectively)
  • Set exchange attempt frequency to 1-2 ps⁻¹

• Production Simulation:

  • Run parallel simulations at assigned temperatures using Langevin dynamics with friction coefficient 0.1-1.0 ps⁻¹
  • Attempt configuration exchanges between adjacent temperature replicas every 100-500 steps using Metropolis criterion: P = min{1, exp[(βi - βj)(Ui - Uj)]}
  • Continue simulation until multiple folding/unfolding events observed at lowest temperature

• Analysis:

  • Calculate free energy surfaces using weighted histogram analysis method (WHAM)
  • Determine folding temperature from specific heat peak
  • Identify transition states from free energy barriers
  • Compute folding rates from transition path sampling

Troubleshooting:

• Poor exchange rates: Adjust temperature distribution to ensure 10-20% acceptance
• Lack of folding events: Increase simulation time or adjust temperature range
• Structural instability: Strengthen prior potentials or reduce time step

Protocol 2: Machine Learning CG Potential with Enhanced Sampling

Objective: Train transferable neural network potential for multiple proteins using enhanced sampling for data generation.

Materials and Setup:

• Reference Data: Unbiased all-atom MD trajectories (if available) or enhanced sampling simulations
• CG Mapping: Backbone-heavy representation (Cα, C, O, N) or Cα-only for larger systems
• NNP Architecture: Deep Graph Network or Behler-Parrinello type neural network
• Training Framework: TensorFlow or PyTorch with molecular dynamics interface

Procedure:

• Data Generation:
  • Perform all-atom enhanced sampling (metadynamics, replica exchange) for 12+ proteins with diverse secondary structures
  • Apply bias along carefully selected collective variables (radius of gyration, native contacts, etc.)
  • Harvest configurations and forces from unbiased regions of trajectory or reweight biased simulations
  • Map all-atom configurations to CG representation using linear mapping operator Ξ

• Network Training:

  • Initialize network parameters with Xavier or He initialization
  • Normalize input features (distances, angles) to zero mean and unit variance
  • Minimize force-matching loss with L2 regularization: L(θ) = (1/3nM) Σc=1M ∥ΞF(rc) + ∇U(Ξrc;θ)∥² + λ∥θ∥²
  • Use Adam optimizer with learning rate decay from 10⁻³ to 10⁻⁶
  • Validate on held-out protein sequences not in training set

• Simulation and Validation:

  • Run CG-MD simulations using trained NNP starting from extended and native states
  • Compute free energy surfaces from accumulated trajectories
  • Compare with experimental data (B-factors, FRET, NMR) when available
  • Assess transferability by simulating mutant sequences

Troubleshooting:

• Poor energy conservation: Reduce time step or increase network stability
• Overfitting: Increase regularization or training set diversity
• Lack of transferability: Include more diverse protein folds in training set

Workflow for ML-CG Potential Development

The Scientist's Toolkit: Research Reagent Solutions

Table 3: Essential Computational Tools for Enhanced Sampling in CG Protein Folding

Tool Category	Specific Implementation	Primary Function	Application Context
CG Force Fields	Martini [6], AWSEM [6], Structure-Based Models [6]	Provide effective interactions between CG sites	Baseline potentials for CG-MD simulations
Neural Network Potentials	Behler-Parrinello NN [6], Deep Graph Networks [53]	Learn many-body interactions from reference data	ML-enhanced CG with quantum accuracy
Enhanced Sampling Packages	PLUMED [52], SSAGES [52]	Implement advanced sampling algorithms	Collective variable-based enhanced sampling
MD Simulation Engines	GROMACS [52], NAMD [52], CHARMM [51], OpenMM	Perform high-performance molecular dynamics	Production simulations with various samplers
Analysis Tools	MDTraj, MDAnalysis, PyEMMA, WHAM	Process trajectories and compute observables	Free energy calculation and state identification
Mapping Optimization	Essential Dynamics CG [7], Mapping Entropy Minimization [53]	Determine optimal CG representations	System-specific coarse-graining

Concluding Remarks

Enhanced sampling techniques have transformed coarse-grained molecular dynamics from a method limited by kinetic traps to a powerful tool for exploring protein folding landscapes. The integration of traditional statistical mechanical approaches with modern machine learning has created a robust framework for addressing the fundamental challenges of rough energy landscapes and rare events. As these methods continue to mature, their capacity to illuminate the relationship between protein sequence, structure, and function will undoubtedly expand, offering unprecedented insights into the molecular mechanisms underlying biological processes and disease states. The protocols and analyses presented here provide researchers with practical strategies to implement these advanced techniques in their own investigations of protein folding dynamics.

The development of machine-learned force fields (ML-FFs) represents a paradigm shift in molecular simulation, offering the potential to bridge the accuracy of quantum mechanics with the computational efficiency of molecular mechanics. However, a significant challenge persists: overfitting. This occurs when a model learns the specific details of its training data—including statistical noise—so closely that it fails to generalize to new, unseen molecular systems or conformations. Within the context of coarse-grained molecular dynamics (MD) for small protein folding, an overfit force field might accurately reproduce the energy landscape of proteins in its training set but fail to predict folding dynamics for novel sequences or stabilize native states correctly, leading to non-physical predictions and unreliable scientific conclusions [55].

The core of the problem lies in the complexity and high-dimensionality of ML-FFs. Models with millions of parameters can easily memorize data without learning the underlying physical principles governing protein interactions. For coarse-grained models, which are often trained on data from all-atom simulations or experimental observables, this can manifest as an incorrect ranking of folded state stability, a bias towards misfolded or non-native states, or an inability to capture the subtle balance of interactions that drive folding [56] [6]. Consequently, implementing robust regularization strategies is not merely a technical exercise but a fundamental requirement for developing transferable, predictive, and reliable ML-FFs for protein science and drug development.

Regularization Strategies for Robust Force Fields

Regularization techniques impose constraints on a machine learning model to prevent it from becoming overly complex. The following strategies are particularly relevant for machine-learned force fields.

Physical Priors and Symmetry Enforcement

Building physical knowledge directly into the model architecture is one of the most powerful forms of regularization. This constrains the model to a physically plausible subspace of functions, drastically reducing its capacity to learn non-physical relationships.

• Energy Conservation and Equivariance: A proper force field must conserve energy and produce forces that are equivariant to rotation and translation (E(3)-equivariance). Models that learn a potential energy surface and derive forces as its negative gradient inherently respect energy conservation. Utilizing E(3)-equivariant neural networks ensures that the model's predictions transform correctly under these operations, a critical physical prior [57] [5].
• Permutation Symmetries: The energy contribution of bonded interactions, such as bonds, angles, and torsions, must be invariant to the permutation of atoms that leave the corresponding internal coordinate unchanged. For instance, the parameters for a bond between atom i and j should be symmetric: ξ(bond)ij = ξ(bond)ji. Enforcing such symmetries in the model architecture reduces the number of free parameters and guarantees physically sensible outputs [57].
• Functional Form Priors: Rather than predicting energies directly from atomic coordinates in a "black box" manner, some ML-FFs learn the parameters for a pre-defined molecular mechanics energy function (e.g., harmonic bonds, periodic torsions). The Grappa model, for example, uses a graph neural network to predict molecular mechanics parameters, thereby restricting the energy function to a known, interpretable form that is computationally efficient and inherently stable [57].

Bayesian Inference and Uncertainty Quantification

Bayesian methods offer a probabilistic framework for regularization by treating model parameters as distributions rather than fixed values. This allows for the explicit quantification of predictive uncertainty.

• BICePs (Bayesian Inference of Conformational Populations): This approach uses Bayesian inference to reweight a simulated conformational ensemble (the prior) against experimental data. It introduces a likelihood function that includes uncertainty parameters (σ) for the experimental data. By sampling the posterior distribution, BICePs simultaneously learns reweighted populations that better agree with experiment and the most likely values for the data uncertainties. The resulting BICePs score serves as a robust metric for force field selection and validation, helping to identify models that generalize well without overfitting to specific data points [58].

Conventional ML Regularization in Loss Functions

Standard machine learning regularization techniques can be incorporated directly into the optimization process.

• L1 (Lasso) and L2 (Ridge) Regularization: These techniques add a penalty term to the loss function during training. L2 regularization penalizes the sum of squared weights (α * Σ||w||²), encouraging the model to keep all weights small but non-zero. L1 regularization penalizes the sum of absolute weights (α * Σ|w|), which can drive many weights to exactly zero, effectively performing feature selection. The hyperparameter α controls the strength of the regularization [59] [60].
• Data Augmentation: While more common in other ML domains, the concept can be applied to ML-FFs by augmenting training data with physically plausible perturbations, such as small rotations or translations, or by training on a diverse set of molecular systems to encourage broader transferability [6].

The table below summarizes the primary regularization strategies and their applications in force field development.

Table 1: Summary of Regularization Strategies for Machine-Learned Force Fields

Strategy	Mechanism	Primary Application in ML-FFs
Physical Priors & Symmetries	Constrains model architecture to obey physical laws (e.g., energy conservation, symmetry).	Ensuring rotational equivariance, permutation invariance of energy terms, and functional form constraints [57] [5].
Bayesian Inference (e.g., BICePs)	Treats parameters as distributions; introduces uncertainty parameters to prevent overconfidence.	Reweighting conformational ensembles and validating force fields against noisy experimental data [58].
L1/L2 Regularization	Adds parameter norm penalties to the loss function to discourage model complexity.	Preventing overfitting in the neural network components of the force field during training [59] [60].
Diverse Training Data	Exposes the model to a wide variety of molecular systems and conformations.	Promoting transferability across protein sequences and folds, preventing memorization of specific systems [6].

The following diagram illustrates how these different regularization strategies are integrated into a typical ML-FF development workflow.

Application Notes: Regularization in Action for Protein Folding

The theoretical strategies outlined above are already being applied to develop next-generation coarse-grained models for protein folding. The following examples showcase their practical implementation and impact.

Case Study: Grappa - A Machine-Learned Molecular Mechanics Force Field

Grappa exemplifies the use of physical priors and symmetry enforcement as a core regularization mechanism. Instead of being a "black box" potential, Grappa employs a graph neural network to predict parameters for a standard molecular mechanics force field (bonds, angles, torsions, etc.). This architectural choice heavily regularizes the model by restricting the energy function to a known, physics-based form. Furthermore, Grappa's "symmetric transformer" explicitly encodes the required permutation symmetries for different interaction types (e.g., ξ(torsion)ijkl = ξ(torsion)lkji), ensuring physical consistency and reducing the risk of learning spurious, non-symmetric correlations [57].

Result: Grappa demonstrates state-of-the-art accuracy on small molecules, peptides, and RNA. It successfully reproduced experimentally measured J-couplings and improved the calculated folding free energy of the small protein chignolin. Its performance shows that a highly regularized model can be both accurate and computationally efficient, operating at the same speed as traditional force fields in MD engines like OpenMM and GROMACS [57].

Case Study: A Universal Coarse-Grained Force Field

This model, described by machine learning coarse-grained potentials of protein dynamics, highlights the importance of diverse training data and thermodynamic consistency for preventing overfitting. The model was trained on a massive, unbiased all-atom MD dataset containing approximately 9 milliseconds of simulation data across twelve different proteins with varied secondary structures. This diversity forces the model to learn generalizable principles of protein energetics rather than memorizing specific proteins [6].

The training used a variational force-matching approach, which aims to make the coarse-grained model thermodynamically consistent with the reference all-atom data. The loss function L(θ) = (1/3nM) Σ || ΞF(r_c) + ∇U(Ξr_c; θ) ||² ensures the CG forces match the mapped all-atom forces. To further regularize the model and prevent unphysical conformations, the authors incorporated simple prior potentials for bonds, repulsive interactions, and chirality. This "delta-learning" approach simplifies the task for the neural network [6].

Result: The resulting single potential was able to simulate all twelve training proteins and, crucially, also captured experimental structural features of mutated proteins not seen during training. This demonstrates successful generalization and a lack of overfitting, enabled by the broad training data and physical priors [6].

Case Study: BICePs for Force Field Selection and Validation

The BICePs method provides a Bayesian framework to regularize and validate force fields after simulation. In one study, conformational ensembles of chignolin generated by nine different force fields were reweighted against 158 experimental NMR measurements. BICePs acts as a regularizer by inferring the most likely uncertainties (σ) in the experimental data and the conformational populations that agree with the data within those uncertainties, preventing overfitting to potentially noisy measurements [58].

Result: The BICePs score provided a clear metric for model selection, identifying which of the nine force fields produced ensembles most consistent with experiment. In all cases, reweighted populations correctly favored the folded state, even for a force field (A99SB-ildn) that initially favored a misfolded state. This shows how Bayesian regularization can be used to assess and correct for force field bias [58].

Table 2: Key Experimental Validation Metrics from Case Studies

Model / Method	Key Validation Metric	Reported Performance	Evidence of Generalization
Grappa [57]	Folding Free Energy (Chignolin)	Improved agreement with experiment	Accurate forces for peptides, RNA, and novel peptide radicals.
Universal CG Model [6]	Structural States & Mutant Energetics	Identified relevant states with comparable energetics	Single model worked for 12 proteins; captured features of mutated proteins.
BICePs [58]	BICePs Score vs. NMR Data	Ranked 9 force fields; identified best agreement	Reweighted ensembles consistently favored native state across force fields.

Experimental Protocols

This section provides detailed methodologies for implementing and validating the regularization strategies discussed.

Protocol: BICePs-Assisted Force Field Validation

This protocol outlines the steps for using BICePs to reweight a conformational ensemble and compute a validation score [58].

• Input Data Preparation:

  • Conformational Ensemble: Generate a Markov State Model (MSM) or a set of conformational states from MD simulations (e.g., performed on Folding@home).
  • Experimental Data: Compile a set of experimental observables and their measured values (e.g., d_j for NMR chemical shifts, J-couplings, or NOE distances). A typical set may include 100+ measurements.

• Forward Model Calculation:

  • For each conformation in the ensemble, calculate the predicted value of each experimental observable (d_j*). For example, use a chemical shift prediction tool like SHIFTX2 for NMR data.

• BICePs Sampling:

  • Define the prior distribution p(X) as the population weights from your MSM or simulation.
  • Define a prior for the uncertainty parameters, p(σ), for each experimental restraint.
  • Use the BICePs algorithm to sample the posterior distribution: p(X, σ | D) ∝ p(D | X, σ) · p(X) p(σ).
  • Run the sampler until convergence to obtain the reweighted populations X and the posterior distributions of the uncertainties σ.

• Analysis and Scoring:

  • Reweighted Ensemble: Analyze the conformational states with the highest posterior populations.
  • BICePs Score: Calculate the BICePs score from the samples. This score represents the free energy cost of applying the experimental restraints and is used for model selection—a better (lower) score indicates a force field whose prior ensemble is more consistent with the experimental data.

Protocol: Training a Regularized ML-FF with Physical Priors

This protocol describes the process for training a machine-learned coarse-grained potential, incorporating physical priors and L2 regularization [6].

• Data Generation and Coarse-Graining:

  • Run extensive all-atom MD simulations for a diverse set of proteins to create a reference dataset of coordinates and forces.
  • Define a coarse-grained mapping (e.g., one bead per Cα atom). Apply this mapping to all atomistic coordinates to get CG coordinates x_c = Ξ r_c.
  • Map the all-atom forces to the CG level: F_CG = Ξ F_AA.

• Model and Loss Function Definition:

  • Architecture: Choose a neural network architecture (e.g., SchNet, Graph Network) to represent the potential energy U(x; θ).
  • Loss Function: Implement the force-matching loss with L2 regularization: L(θ) = (1/3nM) Σ || F_CG + ∇U(x_c; θ) ||² + α ||w||² Here, α is the L2 regularization hyperparameter and w are the network weights.
  • Prior Potentials: Add simple analytical prior terms to the total energy to prevent unphysical states:
    • Bonded Prior: Harmonic potential between consecutive beads to maintain chain connectivity.
    • Repulsive Prior: Repulsive potential (e.g., Lennard-Jones) to prevent bead clashes.
    • Chirality Prior: Dihedral potential to enforce correct chirality and prevent mirror images.

• Model Training:

  • Use a stochastic optimizer (e.g., Adam) to minimize the loss function L(θ) with respect to the network parameters θ.
  • Perform hyperparameter tuning for the regularization strength α and learning rate, using a held-out validation set of proteins.

• Validation and Production Simulation:

  • Equilibrium Properties: Run CG-MD simulations using the trained potential and compute equilibrium properties (e.g., free energy surfaces, folding stability) for proteins not in the training set.
  • Experimental Comparison: Validate simulation outcomes against experimental data, such as folding pathways, native state stability, or NMR observables.

The workflow for this protocol is summarized in the diagram below.

The Scientist's Toolkit

This table lists key software and computational resources essential for implementing the protocols and developing regularized ML-FFs.

Table 3: Essential Research Reagents and Software Solutions

Tool / Resource	Type	Function in Protocol	Key Feature for Regularization
GROMACS [57], OpenMM [57]	MD Software	Running reference all-atom simulations and production CG simulations.	High-performance engines for generating training data and deploying ML-FFs.
BICePs [58]	Software Package	Bayesian reweighting of ensembles against experimental data.	Infers uncertainties and provides a score to prevent overfitting to experimental noise.
Grappa [57]	ML-FF Model	An example of a machine-learned molecular mechanics force field.	Embodies symmetry preservation and functional form priors as core regularization.
PyTorch / TensorFlow	ML Framework	Building and training neural network potentials.	Built-in support for L1/L2 regularization and custom loss functions.
Folding@home [58]	Distributed Computing	Generating massive conformational ensembles for training and validation.	Provides large-scale data to combat overfitting via data diversity.
ESPALOMA [57]	ML-FF Model	A machine-learned force field for organic molecules; a precursor to Grappa.	Demonstrates graph-based parameter assignment.

In the study of protein folding, molecular dynamics (MD) simulations are essential tools for probing thermodynamic and kinetic properties. Traditional simulations often rely on pairwise potentials—models that calculate forces as sums of interactions between two particles. While computationally efficient, these two-body approximations fail to capture critical physical phenomena where the simultaneous interaction between three or more particles generates effects not explainable by pair sums. Multi-body interactions (also called higher-order or N-body interactions) are essential for modeling angular-dependent forces such as electron polarization, hydrogen bonding networks, and cooperative stabilization in protein transition states. In coarse-grained (CG) MD—where groups of atoms are represented as single "beads" to accelerate simulation—capturing these effects is particularly challenging yet critical for maintaining thermodynamic accuracy.

The inclusion of multi-body terms allows CG models to reproduce properties emergent from all-atom simulations, such as free energy landscapes and native contact stability. This document provides application notes and protocols for implementing beyond-pairwise potentials in CG-MD studies of small protein folding, outlining theoretical foundations, practical implementation, and validation techniques tailored for research scientists and drug development professionals.

Theoretical Foundations: Why Move Beyond Pairwise?

The Physical Significance of Multi-Body Forces

In molecular systems, many biologically critical interactions are inherently multi-body. For example:

• Water networks require three-body potentials to reproduce correct structure and dynamics, as the energy of a water molecule in a hydrogen-bonded network depends on the positions of multiple neighbors simultaneously.
• Protein allostery and cooperative folding often involve concerted motions and interactions among three or more residues.
• Electron density redistribution in response to the local environment (polarization) is a fundamental quantum mechanical effect that depends on the configuration of multiple atoms.

Pairwise potentials like Lennard-Jones or Coulombic functions model the potential energy ( \Phi ) as: [ \Phi = \sum{i=1}^{N} \sum{j>i}^{N} \phi{2}(i,j) ] where ( \phi{2}(i,j) ) is the interaction energy between particles i and j. In contrast, a potential including three-body terms becomes: [ \Phi = \sum{i=1}^{N} \sum{j>i}^{N} \phi{2}(i,j) + \sum{i=1}^{N} \sum{j>i}^{N} \sum{k>j>i}^{N} \phi{3}(i,j,k) ] Here, ( \phi{3}(i,j,k) ) represents the non-additive component of the interaction among the triplets of particles i, j, and k. This term can capture angle-dependent forces, which are crucial for maintaining secondary structure elements like alpha-helices and beta-sheets in proteins [61].

The Computational Challenge

The primary challenge with multi-body potentials is their steep computational cost. While pairwise interactions scale as ( \mathcal{O}(N^2) ) with system size N, three-body interactions scale as ( \mathcal{O}(N^3) ), making them prohibitively expensive for large systems [61]. In practice, cutoff distances are applied to limit the number of considered triplets, and advanced algorithms are required to make the calculations tractable for biologically relevant systems like small proteins.

Computational Approaches for Multi-Body Interactions

Several computational strategies have been developed to incorporate multi-body effects into CG-MD simulations of proteins. The table below summarizes the key methodologies, their underlying principles, and applications in protein folding studies.

Table 1: Computational Approaches for Incorporating Multi-Body Interactions in CG-MD

Methodology	Core Principle	Representative Examples	Application in Protein Folding
Explicit N-Body Potentials	Directly parameterize analytical energy functions for particle triplets or quartets.	Axilrod-Teller-Muto (ATM) potential [61].	Rarely used directly for proteins due to high cost; more common in materials science.
Neural Network Potentials (NNPs)	Use machine learning to infer a potential energy function that implicitly includes many-body terms.	Graph Neural Networks (GNNs) trained on all-atom data [62] [6].	Learning a thermodynamically consistent CG potential from all-atom simulations of multiple proteins [6].
Dynamic Elastic Network Models (DENM)	Represent native contacts as breakable springs; unfolding is simulated by breaking non-covalent bonds probabilistically under force.	Extension of standard ENM with a Bell-model based bond breaking probability [7].	Simulating force-induced unfolding of proteins like Fibrinogen and Titin [7].
Explainable AI (XAI) Decomposition	Apply interpretation techniques to "open the black box" of NNPs and quantify learned n-body contributions.	Layer-wise Relevance Propagation for GNNs (GNN-LRP) [62].	Pinpointing stabilizing/destabilizing interactions in metastable protein states and the effects of mutations [62].

Neural Network Potentials for Coarse-Graining

Machine learning, particularly through NNPs, offers a powerful and flexible framework for capturing complex multi-body interactions without needing an explicit functional form.

Theoretical Basis: The core idea is to train a neural network to approximate the Potential of Mean Force (PMF) of the CG system. This is often done via variational force-matching, where the NNP is trained to match the forces derived from all-atom reference simulations [6]. The loss function ( L(\boldsymbol{\theta}) ) for a parameter set ( \boldsymbol{\theta} ) is: [ L(\boldsymbol{\theta}) = \frac{1}{3nM} \sum{c=1}^{M} \parallel \boldsymbol{\Xi} \mathbf{F}(\mathbf{r}c) + \nabla U(\boldsymbol{\Xi} \mathbf{r}c; \boldsymbol{\theta}) \parallel^2 ] where ( \boldsymbol{\Xi} ) is the mapping from all-atom (( \mathbf{r} )) to CG coordinates (( \mathbf{x} )), ( \mathbf{F}(\mathbf{r}c) ) are the all-atom forces, and ( U(\mathbf{x}_c; \boldsymbol{\theta}) ) is the NNP-predicted CG energy [6]. The network architecture (e.g., Graph Neural Networks) inherently learns to represent the multi-body effects present in the training data.

Application Notes: Implementing a Multi-Body CG-MD Workflow

Protocol: Building a Data-Driven Coarse-Grained Potential

This protocol outlines the steps for creating a multi-body-aware CG model for a small protein using neural network potentials, based on the methodology demonstrated in [6].

Step 1: Generate All-Atom Reference Data

• Objective: Create a high-quality dataset of unbiased all-atom MD simulations.
• Procedure:
  • Select a set of small, fast-folding proteins (e.g., Chignolin, WW Domain, Protein G, NTL9) with diverse secondary structures.
  • For each protein, run multiple all-atom MD simulations starting from unfolded or extended states, allowing the proteins to fold and unfold repeatedly.
  • Aggregate simulation data to achieve microsecond- to millisecond-scale cumulative sampling across all proteins. As an example, a recent study accumulated approximately 9 milliseconds of data across 12 proteins [6].
  • Extract and save coordinate-force pairs (( \mathbf{r}c, \mathbf{F}(\mathbf{r}c) )) at regular intervals from the trajectories.

Step 2: Define the Coarse-Grained Mapping

• Objective: Reduce the system's dimensionality by selecting which atoms will be represented by CG beads.
• Procedure:
  • A common and effective mapping is to use a one-bead-per-residue model, where the bead is placed at the Cα atom position.
  • Alternatively, a multi-bead model can be used, representing, for example, the backbone and side chain centers of mass for higher resolution.

Step 3: Design and Train the Neural Network Potential

• Objective: Train an NNP to predict the CG energy and forces from the bead positions.
• Procedure:
  • Choose an Architecture: A Graph Neural Network (GNN) is well-suited as it operates directly on the particle connectivity graph.
  • Implement a "Prior" Potential: Incorporate simple analytical potentials for bonded terms (bonds, angles, dihedrals) and repulsive non-bonded terms to prevent chain rupture and atom clashing. This prior handles basic geometry and lets the NNP learn the complex non-bonded multi-body corrections [6].
  • Training: Minimize the force-matching loss function (shown in Section 3.1) using the all-atom dataset. The NNP learns to predict the forces that would be exerted on the CG beads based on the full atomistic system.

Step 4: Validate the CG Model

• Objective: Ensure the CG model reproduces key thermodynamic and structural properties.
• Procedure:
  • Run CG-MD simulations using the trained NNP.
  • Compare the population of native states and metastable intermediates with all-atom simulations and experimental data (e.g., from NMR).
  • Calculate the free energy landscape as a function of relevant reaction coordinates (e.g., RMSD, radius of gyration). A successful model will show a global minimum at the native state with energetics comparable to the all-atom system [6].
  • Test for transferability by simulating protein mutants (sequences not in the training set) and check if the model predicts stable native-like structures [6].

The following workflow diagram summarizes this multi-step process:

Protocol: Interpreting the NNP with Explainable AI

A key criticism of NNPs is their "black box" nature. The following protocol, based on [62], allows researchers to decompose the model's predictions into physically meaningful n-body contributions.

Step 1: Train a GNN-Based Potential

• Train a CG NNP as described in Protocol 4.1, using a GNN architecture.

Step 2: Apply Layer-wise Relevance Propagation (LRP)

• Objective: Decompose the total predicted energy into contributions from subgraphs (e.g., pairs, triplets) of CG beads.
• Procedure:
  • Use the GNN-LRP technique to attribute the network's energy output to "walks" (sequences of edges) on the input graph [62].
  • Aggregate the relevance scores of all walks associated with a particular subgraph (e.g., a triplet of beads i, j, k) to compute its n-body relevance contribution to the total energy.

Step 3: Analyze Contributions

• The calculated n-body contributions reveal which specific interactions (pairwise and higher-order) the model deems most relevant for stabilizing a given conformation. This can be used to, for example, pinpoint key stabilizing interactions in a folded protein's core or identify which interactions change most significantly upon mutation [62].

The Scientist's Toolkit: Essential Research Reagents & Computational Tools

Successful implementation of multi-body CG-MD relies on a suite of software tools and computational resources. The table below details key components of the research toolkit.

Table 2: Essential Computational Tools for Multi-Body CG-MD Research

Tool Name	Type/Category	Primary Function in Workflow
CHARMM/NAMD	All-Atom MD Engine	Generating reference all-atom simulation data for training. Utilizes force fields like CHARMM36 [63].
RoseTTAFold	AI-based Structure Prediction	Can be used to model missing regions in experimental protein structures prior to simulation [63].
AutoPas	HPC Particle Simulation Library	Provides efficient, parallelized algorithms for calculating two-body and three-body interactions in MD, including novel cutoff methods [61].
GNN-LRP	Explainable AI (XAI) Tool	Decomposes the predictions of a Graph Neural Network potential into human-interpretable n-body contributions [62].
r-RESPA	Multiple Time-Stepping Algorithm	An HPC technique that calculates expensive three-body forces less frequently than two-body forces, improving simulation efficiency [61].

Analysis & Validation: Quantitative Assessment of Multi-Body Effects

Validating a multi-body CG model involves demonstrating its improved accuracy over pairwise models and its fidelity to all-atom and experimental data. Key performance metrics are summarized in the table below.

Table 3: Key Metrics for Validating Multi-Body CG Models in Protein Folding

Validation Metric	Description	Target Outcome	Supporting Data from Literature
Free Energy Landscape	Landscape projected onto order parameters like RMSD.	A clear global minimum at the native state with barrier heights matching reference data.	NNPs recovered experimental conformations and free energy landscapes for 12 tested proteins [6].
Native State Stability	Ability of the model to maintain the folded structure during simulation.	High stability (low RMSD fluctuation) in long simulations.	Including expected persistent interaction features reduced false positive rates in complex classification by 2-5 fold [63].
Transferability	Performance on proteins or mutants not included in training.	Accurate prediction of stable folded structures for new sequences.	A single multi-protein NNP captured experimental structural features of mutated proteins [6].
Interpretability	Ability to attribute energy to specific n-body interactions.	Identification of physically plausible, key stabilizing interactions (e.g., specific salt bridges, hydrophobic clusters).	GNN-LRP interpretation of a CG model for protein NTL9 pinpointed stabilizing/destabilizing interactions in metastable states [62].

The computational study of protein folding is a cornerstone of modern biophysics, with significant implications for understanding cellular function and drug development. Within this field, coarse-grained (CG) molecular dynamics (MD) simulations have emerged as a powerful technique, offering a balance between computational efficiency and biological accuracy. These models simplify atomic detail to focus on essential interactions that govern folding kinetics and thermodynamics. A critical challenge, however, lies in the fact that different protein classes, defined by their dominant secondary structures such as α-helices or β-sheets, present distinct folding landscapes and energetic constraints. This application note details specialized protocols for optimizing CG-MD parameters to accurately simulate the folding of small proteins dominated by either α-helical or β-sheet content, framed within the context of advanced machine-learned and physics-based CG models. We provide a structured comparison of intrinsic conformational preferences, detailed parameterization methodologies, and validated experimental workflows to guide researchers in tailoring their simulations for specific protein structural classes.

Theoretical Foundations and Intrinsic Propensities

Understanding the innate conformational preferences of amino acids is fundamental to parameterizing any protein model. Simplified models that isolate local interactions reveal that steric constraints and backbone stereochemistry create intrinsic biases distinct from frequencies observed in static protein databases.

Hard-Sphere Models and Conformational Equilibrium

Molecular dynamics simulations of an Ala dipeptide mimetic, employing only steric interactions and local stereochemical constraints, demonstrate that β-sheet structures are roughly three and a half times more probable than α-helical structures [64]. This intrinsic preference stands in stark contrast to the overrepresentation of α-helices in the Protein Data Bank, highlighting the influence of non-local interactions in native structures. The equilibrium populations from these hard-sphere simulations show α-helical conformations at less than 25%, while β-sheet and polyproline II (PPII) conformations dominate with approximately 68% probability [64]. These findings underscore the necessity for force fields to accurately capture these underlying preferences as a foundation for modeling specific protein classes.

The Critical Role of the Backbone Bond Angle (τ)

Transitions between α-helix and β-sheet conformations are strongly correlated with the backbone bond angle τ (N–Cα–C). Hard-sphere MD simulations indicate that transitions only occur when τ is greater than 110° [64]. Furthermore, the probability distribution of τ for residues sampling the "bridge" region of φ-ψ space (between canonical α-helix and β-sheet basins) is shifted toward larger angles compared to other regions. This interdependence between dihedral angles and bond angles is a critical, yet often underspecified, parameter in many classical force fields. Modern machine-learned coarse-grained models, trained on diverse all-atom simulation data, are better equipped to capture these complex multi-body correlations, leading to improved performance on proteins containing mixed motifs like β-β-α (BBA) folds [5].

Parameter Optimization Strategies

Quantitative Data for Force Field Parameterization

The table below summarizes key quantitative findings from foundational studies that should inform the parameterization of CG models for specific protein classes.

Table 1: Key Quantitative Parameters for CG-MD Optimization

Parameter / Finding	Description	Relevance to Protein Class	Source
α-helix vs β-sheet intrinsic ratio	β-sheet is ~3.5x more probable than α-helix in steric models.	Basis for ensuring force field balance; underlies helical stability challenge.	[64]
Alanine α′ / β′ equilibrium populations	26% in α-helical (α′) region; 68% in β-sheet/PPII (β′) region.	Reference value for validating amino acid backbone propensity in a CG model.	[64]
Bond angle τ critical value	Transitions between α and β conformations require τ > 110°.	Crucial for modeling conformational dynamics and folding pathways.	[64]
BBA folding challenge	CG models show difficulty stabilizing the native state of BBA, a mixed α/β protein.	Benchmark for testing model transferability across structural classes.	[5]
Machine-learned CG speed gain	Machine-learned CG is orders of magnitude faster than all-atom MD.	Enables sufficient sampling for folding/unfolding transitions in optimization cycles.	[5]

Machine-Learned Coarse-Grained Force Fields

A transformative advance is the development of machine-learned transferable CG models. These models are trained on a large and diverse set of all-atom protein simulations using deep learning methods, such as the variational force-matching approach [5]. The resulting force field, for example, CGSchNet, demonstrates the ability to predict metastable states of folded, unfolded, and intermediate structures for proteins not included in the training set, with low sequence similarity [5]. This "bottom-up" approach ensures that the CG model retains a direct link to the accurate physics of all-atom simulations while achieving the computational efficiency required to sample folding events. Its demonstrated success on small fast-folding proteins like chignolin, TRP-cage, and villin headpiece makes it an ideal platform for further optimization toward specific structural classes [5].

Experimental Protocols

Protocol 1: Optimization of CG Parameters for α-Helical Bundles

This protocol is designed for small, α-helix-rich proteins like the villin headpiece or de novo designed proteins like alpha3D.

Workflow Overview:

Step-by-Step Methodology:

• System Preparation and Training Data Curation:

  • Select a target α-helical protein (e.g., PDB: 1YRF, 2A3D).
  • If using a machine-learned force field like CGSchNet, verify that the target sequence has <40% similarity to the model's training set to ensure genuine predictive power [5].
  • Generate reference data from all-atom simulations if available, focusing on helical stability and native contact formation.

• Model Selection and Prior Energy Term Tuning:

  • Employ a machine-learned CG model known to handle helical structures effectively (e.g., as demonstrated for villin headpiece and alpha3D [5]).
  • For non-machine-learned CG models, strengthen the hydrogen bonding potential and torsional potentials in the φ-ψ regions corresponding to the α-helical basin (approximately φ = -60°, ψ = -45°) to counter the intrinsic under-stabilization suggested by hard-sphere models [64].

• Simulation and Validation:

  • Initiate simulations from extended conformations and monitor the formation of helical bundles.
  • Run parallel-tempering or long constant-temperature Langevin simulations to ensure converged sampling of the equilibrium distribution [5].
  • Key Validation Metrics:
    • Fraction of Native Contacts (Q): The folded state should achieve Q > 0.9 [5].
    • Cα Root-Mean-Square Deviation (RMSD): The dominant cluster should have a low Cα RMSD (< 0.2 nm) from the native structure.
    • Helicity: Monitor the per-residue helicity to ensure it matches experimental data (e.g., from CD spectroscopy) and all-atom reference simulations.
    • Terminal Fluctuations: Compare the Cα root-mean-square fluctuations (RMSF) of the CG folded state with all-atom simulations, expecting similar flexibility patterns, particularly at the termini [5].

Protocol 2: Optimization of CG Parameters for β-Sheet Proteins

This protocol is tailored for proteins rich in β-sheets, which present a challenge due to the need for accurate long-range interactions.

Workflow Overview:

Step-by-Step Methodology:

• System Preparation:

  • Select a target β-sheet protein. Small, well-characterized systems are ideal for initial parameterization.
  • For machine-learned models, the same transferable model used for helical systems can be applied, as its strength is generalizability [5].

• Parameterization of β-Sheet Interactions:

  • The key is to ensure the model can correctly capture the multi-body interactions that stabilize β-sheet formation. Machine-learned force fields inherently learn these complex terms [5].
  • In simpler CG models, special attention must be paid to the backbone bond angle τ. Ensure the model allows for its natural distribution and correlation with dihedral transitions, as a flexible τ > 110° is critical for sampling β-sheet-competent states [64].
  • Explicitly parameterize the directional hydrogen bonding and side-chain packing interactions that define anti-parallel and parallel β-sheets.

• Simulation of Folding and Misfolding Pathways:

  • Start simulations from extended and partially folded conformations.
  • Deliberately screen for a newly identified class of entanglement misfolding, where loops trap other protein sections or native entanglements fail to form. This can be done by analyzing simulation trajectories for non-native contacts and persistent, stable structures that deviate from the native state [65]. All-atom simulations support the existence of these misfolds, which can be stable and evade cellular quality control [65].

• Validation Metrics:

  • Structural Fidelity: Assess the Cα RMSD and fraction of native contacts of the predicted folded state.
  • β-Sheet Content: Quantify the formation of correct β-strands and their proper pairing.
  • Misfolding Propensity: Report the frequency and stability of misfolded entangled states observed in simulations [65].
  • Comparison with Mass Spectrometry: For proteins studied experimentally, compare inferred structural changes from simulation with those obtained from techniques like mass spectrometry, which can probe folding landscapes [65].

Table 2: Key Research Reagent Solutions for CG-MD of Protein Folding

Item / Resource	Function / Description	Application in α-helical vs β-Sheet Studies
Machine-Learned CG Force Field (e.g., CGSchNet)	A transferable coarse-grained force field trained on all-atom MD data with deep learning.	Universal starting point for both protein classes; demonstrates extrapolative power on new sequences [5].
All-Atom MD Force Fields (Amber, CHARMM)	High-resolution force fields for generating reference data and validating CG results.	Used to assess limitations in capturing τ-angle correlations and bridge region conformations [64].
Hard-Sphere Dipeptide Model	A simplified model with only steric interactions and stereochemical constraints.	Used to establish intrinsic β-sheet preference and role of bond angle τ [64].
Parallel-Tempering (Replica-Exchange) MD	An enhanced sampling simulation technique.	Essential for obtaining converged folding/unfolding free energy landscapes for both α-helical and β-sheet proteins [5].
AI2BMD System	An AI-based ab initio biomolecular dynamics system using protein fragmentation.	Provides ab initio accuracy for energy/force calculations for large proteins; useful for generating high-quality training data [66].
Pfam Database	A large collection of protein families and alignments.	Source of diverse protein sequences for pre-training universal representation models [67].

Concluding Remarks

The optimization of coarse-grained MD parameters for specific protein classes is a nuanced process that must bridge fundamental physical principles, such as the intrinsic conformational preferences of the polypeptide backbone, and the complex, emergent interactions that define native folds. The advent of machine-learned, transferable coarse-grained models represents a significant leap forward, providing a robust and efficient platform for simulating a wide range of proteins. By following the protocols outlined herein—which emphasize rigorous validation against experimental data and awareness of class-specific challenges like helical stability and β-sheet entanglement—researchers can enhance the predictive power of their simulations. This approach accelerates in silico protein folding studies and provides deeper insights into the misfolding mechanisms implicated in disease, thereby supporting the broader goals of scientific discovery and therapeutic development.

Molecular dynamics (MD) simulation is a powerful computational technique for studying protein folding, but a fundamental trade-off exists between the atomic-level accuracy of a simulation and its computational cost. All-atom (AA) MD simulations provide high fidelity by explicitly modeling every atom but come at an "extreme computational cost," limiting the time and length scales accessible for study [5]. Coarse-grained (CG) MD addresses this by simplifying the molecular representation, replacing groups of atoms with fewer interacting sites or "beads," thereby dramatically increasing simulation speed and extending the accessible scales [16] [68]. This application note, framed within research on coarse-grained parameters for small protein folding, details the quantitative trade-offs, provides protocols for implementation, and outlines key resources for researchers aiming to leverage CGMD in drug development.

Quantitative Analysis of Efficiency-Accuracy Trade-offs

The choice between simulation methods involves balancing the resolution of information against computational resources. The table below summarizes the core trade-offs between All-Atom, traditional Coarse-Grained, and emerging Machine-Learned Coarse-Grained (ML-CG) models.

Table 1: Performance and Accuracy Trade-offs in Protein Folding Simulations

Simulation Model	Representative Resolution	Computational Speed vs. AA MD	Key Strengths	Key Limitations and Inaccuracies
All-Atom (AA) MD [5] [1]	Every atom explicitly represented.	1x (Baseline)	• High-resolution spatial and temporal data.• Gold standard for accuracy when achievable.	• Extreme computational cost.• Limited to microsecond-millisecond timescales for small proteins.
Traditional Coarse-Grained (CG) Models [16] [68]	1 bead per residue or coarser.	Several orders of magnitude faster.	• Access to larger systems and longer timescales.• Smoother free energy surfaces for faster sampling.	• Difficulty modeling multi-body interactions.• Often limited to specific system types (e.g., MARTINI for membranes) [5].
Machine-Learned Coarse-Grained (ML-CG) Models [5]	1 bead per residue (e.g., Cα atoms).	Several orders of magnitude faster.	• Transferable across sequences.• Accurately predicts metastable states, folding/unfolding, and relative folding free energies [5].	• Requires extensive AA data for training.• Performance may vary on complex motifs like anti-parallel β-sheets [5].

A 2025 study demonstrated the promise of ML-CG models, showing they could successfully predict metastable states of folded and unfolded structures, fluctuations of disordered proteins, and relative folding free energies of mutants, all while being several orders of magnitude faster than a comparable all-atom model [5]. However, the study also noted that while the relative positions of free energy minima were correct, the absolute free energy differences between states did not always exactly match the all-atom reference, illustrating a key area of approximation [5].

Experimental Protocols for CGMD Implementation

Protocol 1: Bottom-Up Development of a Machine-Learned CG Force Field

This protocol outlines the steps for creating a transferable, bottom-up CG force field using deep learning, as exemplified by recent research [5].

• Training Data Generation (AA MD):

  • Run a diverse set of all-atom, explicit-solvent MD simulations of small proteins and peptides. The training set should include various folded structures and di/trimers of peptides to ensure diversity [5].
  • Critical Note: The storage requirements for these reference trajectories can be substantial. It is common practice to save only positional data without force labels to reduce file size, which can constrain later training options [69].

• Model Definition and Training (Force Matching):

  • Coarse-Grained Mapping: Define the mapping from all-atom to CG coordinates. A common scheme is to represent each amino acid with a single bead located at the Cα atom [5] [69].
  • Force Matching Objective: Employ a variational force-matching approach. The learning objective is to minimize the mean-squared error between the instantaneous all-atom forces projected onto the CG coordinates and the forces predicted by the ML-CG model [5] [69].
  • Neural Network Training: Train a neural network (e.g., CGSchNet) on the generated dataset to learn the CG force field. The model learns to represent the effective, multi-body physical interactions between the CG sites [5].

• Model Validation on Unseen Sequences:

  • Test the trained CG model on proteins that were not present in the training or validation datasets and have low sequence similarity (<40%).
  • Perform parallel-tempering and constant-temperature Langevin simulations to obtain converged equilibrium distributions.
  • Validate the model by comparing the resulting free energy landscapes, folded state geometries (using metrics like Cα root-mean-square deviation and fraction of native contacts), and protein fluctuations against reference all-atom simulations or experimental data [5].

Protocol 2: Refining General-Purpose CG Topologies with Bayesian Optimization

This protocol is ideal for improving the accuracy of established, general-purpose CG force fields (e.g., Martini3) for specific protein classes, using an active learning framework [37].

• Define Target and Parameter Space:

  • Select target macroscopic properties for optimization, such as protein density (ρ) or radius of gyration (Rg). These are calculated from reference all-atom simulations or obtained from experiments [37].
  • Identify the low-dimensional set of bonded parameters, (\boldsymbol{\theta}), within the CG topology to be optimized. This typically includes bond lengths ((b0)), bond constants ((kb)), angle magnitudes ((\Phi)), and angle constants ((k_\Phi)) [37].

• Set Up Bayesian Optimization (BO) Loop:

  • Initialize Surrogate Model: The BO algorithm uses a probabilistic surrogate model to approximate the relationship between the CG parameters ((\boldsymbol{\theta})) and the target properties [37].
  • Define Acquisition Function: This function, guided by the surrogate model, balances exploration of new parameter regions and exploitation of known promising areas [37].
  • Iterate Until Convergence: a. The acquisition function proposes a new set of parameters to evaluate. b. Run a CG simulation with the proposed parameters. c. Calculate the target properties (ρ, Rg) from the simulation output. d. Update the surrogate model with the new result. e. Repeat until the difference between CG and reference properties is minimized.

• Validation Across Scales:

  • The final optimized parameter set should be validated by simulating systems of varying complexity (e.g., different protein sizes or degrees of polymerization) not used during optimization to ensure transferability [37].

Visualization of Workflows and Logical Relationships

CG Force Field Development Pathways

The diagram below outlines the two primary methodologies for developing coarse-grained models.

The Accuracy vs. Speed Trade-off in Simulation Models

This diagram conceptualizes the fundamental relationship between simulation resolution and computational efficiency.

The Scientist's Toolkit: Research Reagent Solutions

Table 2: Essential Resources for Coarse-Grained Protein Folding Studies

Tool Name	Type	Primary Function	Relevance to CG Protein Folding
CGSchNet [5]	Machine-Learned Force Field	A neural network-based CG model trained via force matching.	Provides a transferable, bottom-up force field for predicting protein folding landscapes and dynamics on new sequences.
Martini 3 [37]	Top-Down Coarse-Grained Force Field	A general-purpose force field parameterized against experimental thermodynamic data.	A widely used model; its bonded topologies can be refined for specific protein classes using optimization frameworks.
CALVADOS [70]	Coarse-Grained Model	A residue-level model optimized for simulating intrinsically disordered proteins (IDPs).	Specialized for studying the folding, dynamics, and phase separation of IDPs, which have highly heterogeneous conformational ensembles.
Bayesian Optimization (BO) Framework [37]	Optimization Algorithm	Refines bonded parameters in CG topologies against target properties (e.g., from AA MD or experiment).	Bridges the accuracy gap of general-purpose models like Martini for specific applications, ensuring both efficiency and accuracy.
Variational Force-Matching [5] [69]	Training Methodology	Learns CG force fields by matching projected all-atom forces.	The foundational principle for many modern, data-driven bottom-up CG force fields, enabling high accuracy.

Molecular dynamics (MD) simulations provide atomic-level insights into protein folding mechanisms, but achieving physically realistic representations remains challenging, particularly for coarse-grained (CG) models that sacrifice atomic detail for computational efficiency [71]. The core challenge lies in ensuring these simplified models avoid unphysical states—conformations with atomic overlaps, distorted geometries, or unrealistic interactions—while preserving accurate folding thermodynamics and kinetics [6]. As CG models increasingly incorporate machine learning potentials trained on all-atom simulation data [6], robust validation protocols become essential for distinguishing physically realistic predictions from computational artifacts. This application note establishes standardized protocols for validating physical realism in CG-MD parameters, specifically focusing on small protein folding studies where accurate energy landscapes are critical for meaningful biological interpretation.

Quantitative Validation Metrics and Criteria

Key Physical Realism Metrics for Protein Folding Simulations

Validation Category	Specific Metric	Target Value	Measurement Method
Steric Clashes	Atomic overlap violations	< 1% of residue pairs [72]	Distance < sum of van der Waals radii [72]
Backbone Geometry	Ramachandran outliers	< 0.05% [72]	φ/ψ dihedral distribution analysis
Bond Geometry	Bond length/angle deviations	< 4× standard deviation [72]	Comparison to Engh-Huber parameters [72]
Energy Landscape	Folding free energy (ΔG)	~0.8 kcal/mol [73]	Dual-cutoff state population analysis [73]
Structural Accuracy	Cα-RMSD to native state	< 1.5 Å [73]	Clustering with 1.0 Å cutoff [73]
Force Field Balance	Folded vs. disordered stability	Maintains both [74]	SAXS/NMR comparison [74]

Experimental Method	Validated Property	Application in CG-MD Validation
X-ray Crystallography	Native-state structure, B-factors	B-factor correlation > 0.7 [7], RMSD validation [72]
NMR Spectroscopy	Secondary structure propensities, dynamics	Chemical shift comparison [74]
Small-Angle X-Ray Scattering (SAXS)	Chain dimensions, ensemble properties	Radius of gyration (Rg) validation [74]
Folding Kinetics	Folding rates, pathways	Folding time comparison [73]
Atomic Force Microscopy (AFM)	Mechanical unfolding pathways	Force-extension curve matching [7]

Experimental Protocols for Physical Realism Validation

Protocol 1: Steric Clash Detection and Resolution

Purpose: To identify and quantify unphysical atomic overlaps in CG-MD simulation trajectories.

Procedure:

• Trajectory Sampling: Extract frames from production MD simulations at regular intervals (e.g., every 100ps for microsecond trajectories) [73]
• Distance Calculation: Compute pairwise atomic distances for all non-bonded atoms within a 5Å cutoff
• Clash Identification: Flag atom pairs where distance < (vdWradiusi + vdWradiusj) × 0.8 [72]
• Quantification: Calculate clash frequency as percentage of residues involved in steric violations per frame [72]
• Spatial Mapping: Identify persistent clash regions in protein structure to pinpoint problematic force field terms

Acceptance Criteria: Successful refinement should reduce steric clashes to <1% of residue pairs while maintaining backbone topology (TM-score >0.8) [72].

Protocol 2: Backbone Geometry Validation

Purpose: To ensure protein backbone conformations sample physically realistic φ/ψ dihedral space.

Procedure:

• Dihedral Calculation: Compute φ/ψ angles for all non-proline/glycine residues across simulation trajectory
• Ramachandran Assessment: Categorize each (φ,ψ) pair as favored, allowed, or outlier using Richardson's criteria [72]
• Temporal Analysis: Track dihedral transitions to identify unstable regions in force field
• Secondary Structure Validation: Compare α-helix and β-sheet propensities to experimental measurements (e.g., NMR chemical shifts) [74]

Troubleshooting: If outlier rates exceed 0.05%, consider adding backbone torsion priors to CG potential as in Wang et al. [6]

Protocol 3: Force Field Transferability Testing

Purpose: To validate CG-MD parameters across multiple protein systems and structural classes.

Procedure:

• Diverse Dataset Selection: Include proteins with varied secondary structure (α-helical, β-sheet, mixed αβ) as in Wang et al.'s 12-protein dataset [6]
• Folding Thermodynamics: Calculate folding free energies for multiple proteins using same parameter set
• Mutational Sensitivity: Test on protein mutants not included in training data [6]
• Comparative Analysis: Verify consistent performance across structural classes and chain lengths

Validation Benchmark: Successful transferability is demonstrated when a single CG potential reproduces experimental structural features across multiple proteins [6].

Workflow Visualization for Validation Protocols

Physical Realism Validation Workflow

Steric Clash Identification Pipeline

Research Reagent Solutions for CG-MD Validation

Essential Software Tools and Datasets

Tool/Dataset	Application	Key Features
ModRefiner [72]	Atomic model refinement from Cα traces	Physics- and knowledge-based force field; improves steric clashes and H-bond networks
Dynamic Elastic Network Model (DENM) [7]	Coarse-grained unfolding simulation	Native structure-based potentials; Bell model for bond breaking
Machine Learning Coarse-Grained Potentials [6]	Neural network CG force fields	Trained on multi-ms all-atom MD; preserves thermodynamics
AMBER/CHARMM Force Fields [73] [74]	All-atom reference simulations	Balanced protein-water interactions; validated folding properties
Multi-Protein MD Dataset [6]	CG potential training/validation	9ms all-atom simulations for 12 proteins with diverse structures
SAXS/NMR Reference Data [74]	Experimental validation	Chain dimensions and secondary structure propensities

Validating physical realism in coarse-grained MD parameters requires a multi-faceted approach addressing steric clashes, backbone geometry, energy landscapes, and experimental consistency. The protocols outlined provide a standardized framework for researchers to identify and rectify unphysical states in protein folding simulations. Special attention should be paid to force field transferability across different protein structural classes and the integration of experimental data throughout the validation process. As machine learning approaches continue to enhance CG-MD methodologies [6], these validation standards will ensure that computational gains in efficiency do not come at the cost of physical accuracy, particularly for critical applications in drug development and protein engineering where reliable atomic-level insights are essential.

Benchmarking CG Performance: Validation Against Experiments and All-Atom MD

In the study of small protein folding using coarse-grained (CG) molecular dynamics (MD), selecting appropriate quantitative metrics is fundamental for accurately characterizing folding pathways, transition states, and thermodynamics. The choice of coordinate directly influences the interpretation of the folding mechanism and the assessment of a model's quality. This application note details three cornerstone metrics—Free Energy Landscapes (FELs), Root-Mean-Square Deviation (RMSD), and Native Contacts (Q)—within the context of parameterizing and validating CG models. We provide a structured comparison of their applications, detailed protocols for their calculation, and visualization of their interplay to equip researchers with a practical toolkit for robust analysis of protein folding simulations.

Quantitative Metrics Comparison

The table below summarizes the core characteristics, strengths, and limitations of the three primary metrics used in analyzing protein folding simulations.

Table 1: Comparative Analysis of Key Quantitative Metrics for Protein Folding Studies

Metric	Definition & Purpose	Key Strengths	Inherent Limitations	Common CG Model Applications
Free Energy Landscape (FEL)	A hypersurface plotting free energy (G) against collective variables (e.g., Q, RMSD). Reveals stable states, barriers, and folding pathways. [75] [76]	Direct visualization of thermodynamics and kinetics; identifies metastable intermediates and transition states. [76]	Computationally expensive to converge; choice of reaction coordinate can distort perception of the landscape. [75]	Essential for validating the thermodynamic consistency of bottom-up CG models. [27] [77]
Root-Mean-Square Deviation (RMSD)	Measures the average displacement of atomic coordinates relative to a reference structure, typically the native state. [76]	Simple, intuitive measure of global structural similarity; widely used and easy to compute.	Global measure insensitive to localized structural changes; results can be skewed by flexible termini or loops. [78]	Standard metric for assessing structural fidelity and fold stability in simulations. [78] [74]
Fraction of Native Contacts (Q)	The fraction of residue-residue contacts present in the simulation that are also present in the native structure. [79]	Excellent reaction coordinate for folding; strongly correlates with folding mechanism in both all-atom and Gō models. [79]	Definition of a "contact" can vary; inherently blind to non-native interactions.	Foundation of structure-based Gō models; primary metric for defining folding progress. [75] [79]

Experimental Protocols & Workflows

Protocol: Constructing a Free Energy Landscape

This protocol describes the calculation of a free energy landscape from a molecular dynamics trajectory using collective variables.

1. Prerequisites

• Input Data: A molecular dynamics trajectory (CG or all-atom) that samples the folding and unfolding events.
• Collective Variables (CVs): Select one or two CVs that describe the folding process. The fraction of native contacts (Q) is often an excellent choice. [79]

2. Calculation Steps

• Step 1: For each frame in the trajectory, compute the values of the chosen CVs (e.g., Q and RMSD).
• Step 2: Construct a 2D histogram of the CVs. This represents the probability distribution, P(Q, RMSD).
• Step 3: Calculate the free energy using the relationship: G(Q, RMSD) = -kₚT * ln(P(Q, RMSD)), where kₚ is Boltzmann's constant and T is the simulation temperature.
• Step 4: Plot the free energy as a contour or surface plot. The minima correspond to stable states (e.g., folded, unfolded), and the saddle points correspond to transition states. [76]

3. Analysis

• Identify the number and location of free energy minima.
• Trace the most probable pathway between minima to infer the folding mechanism.
• The difference in free energy between the folded and unfolded minima yields the folding free energy, ΔG.

Protocol: Calculating Native Contacts (Q)

The fraction of native contacts is a crucial metric for defining the progress of folding.

1. Define Native State Contacts

• Use a representative native structure (e.g., from PDB).
• For each pair of residues (i, j), calculate the distance between heavy atoms or Cα atoms.
• Define a native contact if the distance is below a cutoff (typically 4.5-8.0 Å) and |i-j| > 3 to exclude sequential neighbors. [79]

2. Compute Q During Simulation

• For each simulation frame, calculate the current distance for every predefined native contact pair.
• A contact is considered formed if the current distance is within a specified tolerance (e.g., 1.2-1.5 times) of the native distance. [75] [79]
• Q = (Number of formed native contacts) / (Total number of native contacts).

3. Advanced Application: Transition Path Analysis

• To identify contacts critical for the folding mechanism, calculate the lifetime of each native contact during folding transition paths versus the unfolded state. [79]
• Contacts with a high log-ratio of lifetimes are deemed crucial for the folding mechanism.

Protocol: Utilizing RMSD for Foldability and Stability

RMSD provides a straightforward measure of global structural change.

1. Calculation

• Alignment: Superimpose the protein structure from a simulation frame onto the reference native structure, typically using the Cα atoms of structured regions.
• Calculation: Compute the RMSD as the square root of the average squared distance between the corresponding atoms after alignment.

2. Application in Foldability Assessment

• Recent studies on protein mutants and de novo designs show that the RMSD of predicted structures (e.g., from AlphaFold2) to a reference is quantitatively correlated with protein foldability. [78]
• Higher RMSD values in computational models indicate lower foldability.

The Scientist's Toolkit: Research Reagent Solutions

Table 2: Essential Software and Force Fields for CG Protein Folding Studies

Category	Tool/Reagent	Primary Function	Application Note
CG Simulation Packages	almost (Open-Almost)	Monte Carlo simulation package for all-atom and Cα Gō models. [75]	Enables reversible folding simulations of peptides with various energy functions.
	CHARMM	Molecular dynamics program with CG capabilities. [75]	Used with specific parameter files (e.g., for Cα-KB-GEF model); supports Langevin dynamics.
Force Fields & Models	Gō-like Energy Functions (GEFs)	Structure-based potentials where only native contacts are attractive. [75] [79]	Ideal for studying folding mechanisms due to minimal frustration; available in all-atom and Cα-only variants.
	CGnets (Deep Learning)	Neural networks to learn CG free energy functions from all-atom data. [27]	Captures complex multibody terms in the PMF, overcoming limitations of classical pairwise CG models.
	FLARE (SGP Framework)	Bayesian active learning for on-the-fly CG model development. [77]	Uses sparse Gaussian processes to quantify uncertainty and autonomously collect training data.
Analysis Software	MMTSB Tool Set	Web server and scripts to generate topology/parameters for CG models like Cα-KB-GEF. [75]	Facilitates the setup of sophisticated structure-based simulations.

Logical Workflow and Pathway Visualization

The following diagram illustrates the strategic relationship between the three metrics and a recommended workflow for their application in a CG-MD study of protein folding.

The workflow emphasizes that analysis is not a linear end-point but an integral part of an iterative cycle for model refinement. Validation against experimental data (e.g., folding rates, φ-values) is crucial for establishing the predictive power of the CG model and simulation approach. [79] [76]

The development of coarse-grained (CG) molecular dynamics (MD) models for protein folding represents a significant pursuit in computational biophysics, aiming to bridge the gap between the high accuracy but computational expense of all-atom MD and the need to simulate biologically relevant timescales [5] [80]. A central challenge in this field is achieving true transferability—the ability of a model parameterized on a specific set of proteins to make accurate predictions for novel protein sequences entirely absent from its training data [5]. This application note details the experimental and computational protocols for assessing the transferability of machine-learned coarse-grained potentials, a critical validation step within broader research on CG-MD parameters for small protein folding studies.

Recent advances have demonstrated that machine learning (ML), particularly deep neural networks, can learn transferable CG force fields from large, diverse datasets of all-atom simulations [5] [6]. The table below summarizes the performance of such models on various protein classes unseen during training.

Table 1: Performance of Machine-Learned Coarse-Grained Models on Unseen Proteins

Protein Class/Example	Key Metric	Model Performance on Unseen Sequences	Comparison to Reference
Small Fast-Folding Proteins (e.g., Chignolin, TRPcage, Villin Headpiece)	Folding Free Energy Landscape, Fraction of Native Contacts (Q)	Predicts metastable folded, unfolded, and intermediate states; correctly identifies native state as a low free-energy basin [5].	Closely matches the landscape from all-atom MD reference simulations; comparable population of states [5].
Proteins with Mixed Motifs (e.g., BBA)	Fraction of Native Contacts (Q), Cα Root-Mean-Square Deviation (RMSD)	Predicts native state as a local free-energy minimum, though not always the global minimum, indicating slight inaccuracies in relative stability [5].	Shows greater difficulty in perfectly matching reference landscapes compared to proteins with simpler motifs [5].
Larger, More Complex Proteins (e.g., Engrailed Homeodomain - 1ENH, alpha3D - 2A3D)	Cα Root-Mean-Square Fluctuation (RMSF), Cα RMSD	Folds from extended configurations to native-like structures; stabilizes folded state with flexibility patterns similar to atomistic simulations [5].	Folded state stability and terminal flexibility are consistent with all-atom MD; some sequence regions show higher fluctuations [5].
Intrinsically Disordered Proteins (IDPs)	Conformational Ensemble	Captures the broad conformational distribution of disordered peptides, consistent with NMR data [5] [81].	Reproduces experimental observables better than models lacking transferability [5].
Protein Mutants	Relative Folding Free Energy	Accurately predicts changes in folding free energy due to mutations for sequences outside the training set [6].	Predictions are consistent with experimental measurements, even where converged all-atom simulations are unavailable [5].

Experimental Protocols for Transferability Testing

A robust transferability assessment requires a structured pipeline, from dataset curation to final validation. The following protocol outlines the key stages.

Training Set Curation and Model Learning

The foundation of a transferable model is a diverse and representative training dataset.

• Objective: To assemble a set of all-atom MD simulation data that encompasses a wide variety of protein sequences, secondary structures, and topologies.
• Procedure:
  • Select Training Proteins: Choose proteins with diverse secondary structure arrangements (α-helical, β-sheet, mixed αβ) and varying lengths [6]. The dataset should include monomers and small oligomers (e.g., dipeptides) to sample intermolecular interactions [5].
  • Generate All-Atom Reference Data: Perform unbiased, explicit-solvent all-atom MD simulations for each training protein. The cumulative simulation time should be sufficient to observe folding/unfolding events and sample conformational diversity (e.g., multi-millisecond totals) [6].
  • Define CG Mapping: Reduce the atomic resolution of the training data to a CG representation. A common and effective mapping is to a Cα-only model, where each amino acid is represented by a single bead at the Cα position [5] [7].
  • Train the ML-CG Force Field: Employ a deep learning architecture (e.g., Graph Neural Networks like CGSchNet) and a bottom-up learning objective such as variational force-matching [5] [6].
    • Force-Matching Loss Function: The neural network potential ( U(\mathbf{x}; \boldsymbol{\theta}) ) is trained by minimizing the loss: ( L( \boldsymbol{\theta} ) = \frac{1}{M} \sum{c=1}^{M} \| \boldsymbol{\Xi} \mathbf{F}{AA}(\mathbf{r}c) + \nabla U( \boldsymbol{\Xi} \mathbf{r}c; \boldsymbol{\theta} ) \|^2 ) where ( \boldsymbol{\Xi} \mathbf{F}_{AA} ) are the coarse-grained forces from all-atom data, and ( \nabla U ) are the forces predicted by the neural network potential [6].
  • Incorporate Physical Priors: Augment the machine-learned energy with analytical prior energy terms. These priors enforce basic physics and prevent unphysical states, typically including [5] [6]:
    • Bonded potentials (to maintain chain connectivity)
    • Repulsive potentials (to prevent bead clashes)
    • Chirality restraints (to prevent mirror-image structures)

Rigorous Testing on Unseen Sequences

The core of transferability testing involves evaluating the model on a held-out test set of proteins.

• Objective: To validate that the model can simulate the equilibrium thermodynamics and dynamics of proteins with low sequence similarity to the training data.
• Procedure:
  • Construct a Test Set: Select proteins that have <40% sequence similarity to any protein in the training or validation sets [5]. This ensures a rigorous test of extrapolation rather than interpolation.
  • Perform CG-MD Simulations: Run molecular dynamics simulations using the trained ML-CG force field. To ensure proper sampling of the free energy landscape, especially for folding transitions, use enhanced sampling techniques like Parallel Tempering (Replica Exchange) MD [5].
  • Define Reaction Coordinates: Identify collective variables that describe the folding process. Standard metrics include:
    • Fraction of Native Contacts (Q): Measures the fraction of native pairwise contacts that are formed in a given conformation.
    • Cα Root-Mean-Square Deviation (RMSD): Measures the structural deviation from the known native state.

Validation and Benchmarking

The final step is to quantitatively compare the model's predictions against ground truth data.

• Objective: To benchmark the accuracy and predictive power of the CG model.
• Procedure:
  • Compare to All-Atom MD: For smaller test proteins where feasible, compute the free energy surface as a function of the reaction coordinates (e.g., FES vs. Q and RMSD) from both the CG simulations and reference all-atom MD simulations. Visually and quantitatively assess the overlap of metastable states and the relative free energies of folded, intermediate, and unfolded basins [5].
  • Compare to Experimental Data: For larger proteins where all-atom reference data is unavailable, validate against experimental data. This can include:
    • Cα Root-Mean-Square Fluctuations (RMSF): Compare the flexibility of residues in the CG-folded state to B-factors from X-ray crystallography or fluctuations from all-atom MD [5].
    • Native Structure Stability: Verify that the model's lowest free-energy minimum corresponds to a structure closely resembling the experimental native state (e.g., from PDB) [5].
    • Mutation Effects: Predict the relative folding free energies of mutated proteins and compare the results to experimental measurements [5] [6].

Workflow for Rigorous Transferability Testing of CG Models

The following table catalogues critical computational tools and data resources employed in the development and testing of transferable ML-CG models.

Table 2: Key Research Reagent Solutions for ML-CG Development

Resource Name / Type	Function / Description	Relevance to Transferability Testing
All-Atom MD Datasets (e.g., from [6])	Large-scale, unbiased simulation data of folding proteins used for training.	Provides the foundational "ground truth" for learning a transferable potential; diversity is key to success [5] [6].
CG Mapping Software	Scripts/algorithms to reduce all-atom trajectories to CG representations (e.g., Cα atoms).	Defines the resolution of the CG model and is the first step in the force-matching pipeline [7].
Machine Learning Frameworks (e.g., PyTorch, TensorFlow)	Platforms for implementing and training neural network potentials (NNPs).	Enables the creation of complex, multi-body CG force fields like CGSchNet [5] [6].
Variational Force-Matching Code	Software that implements the force-matching loss function for training.	The core algorithm for ensuring the CG model is thermodynamically consistent with the all-atom data [6].
MD Simulation Engines (e.g., GROMACS, OpenMM)	High-performance software to run CG-MD simulations.	Used to perform production simulations of unseen test proteins to assess model performance [5] [82].
Enhanced Sampling Tools (e.g., for Parallel Tempering)	Modules/plugins for accelerated sampling of free energy landscapes.	Crucial for obtaining converged folding/unfolding statistics for the test proteins during validation [5].
Analysis Suites (e.g., MDTraj, MDAnalysis)	Python libraries for calculating Q, RMSD, RMSF, and other metrics.	Used to quantify and compare the simulation outcomes against references during the benchmarking phase [5].

The protocols outlined herein provide a framework for rigorously evaluating the transferability of machine-learned coarse-grained models. Current evidence demonstrates that ML-based CG force fields, trained on diverse datasets with bottom-up methods, can successfully extrapolate to predict the folding landscapes and thermodynamic properties of proteins unseen during training [5] [6]. This capability, which marks a significant step toward universal CG protein models, enables researchers to explore protein dynamics and folding at a fraction of the computational cost of all-atom simulations, opening new avenues in protein engineering and drug development.

Molecular dynamics (MD) simulation serves as a crucial computational technique for elucidating the physical principles governing protein structures, functions, and dynamics. The most universally predictive protein simulation models traditionally employ all-atom molecular dynamics (AA-MD), which explicitly represents every atom in the system but comes at extreme computational cost [5]. For studying complex processes like protein folding and exploring metastable states, this computational burden often limits the timescales and system sizes that can be practically simulated.

Coarse-grained molecular dynamics (CG-MD) addresses this challenge by representing groups of atoms with single interaction sites, significantly reducing computational complexity. While various CG models have been developed, many have been limited to specific system types or lacked transferability across different protein sequences. Recent advancements combine deep learning methods with diverse training sets from all-atom simulations to develop transferable CG force fields that maintain accuracy while achieving speed increases of several orders of magnitude over atomistic models [5]. This application note examines the capabilities of these next-generation coarse-grained models in capturing metastable states and folding mechanisms compared to traditional all-atom MD, providing practical guidance for researchers investigating protein dynamics and folding.

Performance Comparison: All-Atom vs. Coarse-Grained MD

Quantitative Assessment of Simulation Performance

Table 1: Comparison of All-Atom and Coarse-Grained MD Approaches for Protein Folding Studies

Aspect	All-Atom MD (AA-MD)	Machine-Learned Coarse-Grained MD (CG-MD)
Computational Speed	Reference baseline	Several orders of magnitude faster [5]
Metastable State Prediction	Accurately captures folded, unfolded, and intermediate states [5]	Successfully predicts metastable states comparable to AA-MD [5]
Folding Free Energy Landscapes	Provides reference landscapes [5]	Closely matches AA-MD landscapes for small proteins and peptides [5]
Intrinsically Disordered Proteins	Can capture fluctuations but requires extensive sampling [5]	Accurately models fluctuations and disordered states [5]
Transferability	Universally applicable but computationally expensive [5]	Chemically transferable to new sequences not in training set [5]
Relative Folding Free Energies	Computationally demanding for protein mutants [5]	Predicts relative folding free energies of mutants efficiently [5]
System Size Limitations	Limited to smaller proteins or shorter timescales [5]	Enables simulation of larger proteins (e.g., 54-73 residues) [5]

Application-Specific Performance Insights

For small fast-folding proteins such as chignolin, TRPcage, beta-beta-alpha fold (BBA), and villin headpiece, machine-learned CG models demonstrate remarkable accuracy in predicting metastable folding states. These models stabilize folded states with fraction of native contacts (Q) close to 1 and low Cα root-mean-square deviation values, producing structures that closely resemble the correct native state [5]. The CG model also captures misfolded states observed in reference atomistic simulations, such as the misaligned TYR1 and TYR2 residues in chignolin [5].

For intrinsically disordered proteins (IDPs) and larger protein systems where all-atom simulations struggle to achieve sufficient sampling, CG models demonstrate particular value. They successfully model the fluctuations of disordered regions and enable folding studies of proteins like the 54-residue engrailed homeodomain and 73-residue de novo designed protein alpha3D, which are challenging for conventional AA-MD [5]. When comparing Cα root-mean-square fluctuations within folded states, CG models show similar terminal flexibility to all-atom simulations, though with slightly higher fluctuations along the sequence [5].

Table 2: CG-MD Performance Across Different Protein Classes

Protein Class	Example Systems	CG-MD Performance	Key Metrics
Small Fast-Folding Proteins	Chignolin (2RVD), TRPcage (2JOF)	Accurate prediction of folded states and folding/unfolding transitions [5]	Fraction of native contacts (Q ~1), low RMSD to native state [5]
Mixed Motif Proteins	BBA (1FME)	Captures local minimum but challenges in relative free energy differences [5]	Correct identification of native state as local minimum [5]
Larger Structured Proteins	Engrailed homeodomain (1ENH), alpha3D (2A3D)	Successful folding from extended configurations [5]	Similar terminal flexibility to AA-MD, higher sequence fluctuations [5]
Intrinsically Disordered Proteins	Various disordered peptides	Accurate modeling of fluctuations [5]	Consistent with experimental and AA-MD observations [5]
Protein Mutants	Various folding free energy mutants	Prediction of relative folding free energies [5]	Agreement with experimental trends where AA-MD is unavailable [5]

Experimental Protocols

Workflow for Comparative Folding Studies

The following diagram illustrates the integrated workflow for comparing all-atom and coarse-grained molecular dynamics approaches in protein folding studies:

Protocol 1: All-Atom MD for Folding Studies

System Setup

• Initial Structure Preparation: Obtain starting coordinates from experimental structures (e.g., Protein Data Bank) or generate extended conformations for folding studies [5].
• Force Field Selection: Employ state-of-the-art force fields such as Amber ff19SB [22] for proteins, combined with appropriate water models like OPC [22]. Older force fields (ff99SB, ff14SB) may show greater deviations from experimental data for IDPs [22].
• Solvation and Ionization: Solvate the system in explicit water (e.g., TIP3P water model [83]) using a cubic box large enough to avoid self-interaction through periodic boundaries (e.g., 89 Å × 89 Å × 89 Å for small systems) [83]. Add ions to neutralize system charge.
• Energy Minimization: Perform steepest descent or conjugate gradient minimization to remove steric clashes and bad contacts.

Production Simulation

• Equilibration Protocol: Gradually heat the system to target temperature (typically 300K) using Langevin dynamics or Nosé-Hoover thermostats [83], with gradual release of position restraints on the protein.
• Sampling Strategy: For adequate conformational sampling of folding landscapes, employ enhanced sampling techniques such as replica exchange MD [22] or run multiple independent trajectories [22]. Simulation timescales typically range from nanoseconds to milliseconds depending on system size and folding timescales.
• Parameter Settings: Use a time step of 2 fs [83], with constraints on bonds involving hydrogen atoms (SHAKE algorithm [83]). Maintain constant temperature (303K) and pressure (1 atm) using appropriate thermostats and barostats [83]. Treat long-range electrostatics with particle mesh Ewald summation [83].

Analysis Methods

• Free Energy Landscapes: Construct free energy surfaces using collective variables such as fraction of native contacts (Q) and Cα root-mean-square deviation (RMSD) from native structure [5].
• Metastable State Identification: Use clustering algorithms or Markov state models to identify and characterize metastable states (folded, unfolded, intermediate) [84].
• Validation Metrics: Compare simulation results with experimental data including NMR restraints [84], chemical shifts [22], and kinetic measurements.

Protocol 2: Machine-Learned Coarse-Grained MD

Model Training and Setup

• Training Data Generation: Create a diverse dataset of all-atom explicit solvent simulations of small proteins with various folded structures, along with dimers of mono- and dipeptides [5]. Ensure training covers diverse structural motifs and sequence patterns.
• Neural Network Training: Employ deep learning architectures such as CGSchNet [5] trained using variational force-matching approach [5] to learn effective multi-body interactions from all-atom reference data.
• Transferability Validation: Validate model on proteins with low sequence similarity (16-40%) to training set sequences [5] to ensure generalizability.

Production Simulation

• System Initialization: Initialize CG systems from all-atom structures or extended conformations using mapping functions that reduce atomic resolution to CG beads.
• Sampling Enhancement: Leverage accelerated dynamics from smoothed energy landscapes to sample folding/unfolding transitions. Perform parallel-tempering simulations to ensure converged sampling of equilibrium distribution [5].
• Efficient Propagation: Utilize longer time steps (enabled by reduced resolution) and faster force calculations to achieve microsecond to millisecond effective sampling times.

Analysis and Backmapping

• CG Trajectory Analysis: Analyze simulations using CG-native metrics (e.g., CG distance maps, native contact fractions). Compare folding mechanisms and metastable states with all-atom references [5].
• Backmapping: Where atomistic details are needed, employ reverse coarse-graining (rCG) techniques to reconstruct all-atom coordinates from CG simulations [83] for higher-resolution analysis.
• Experimental Validation: Validate predictions against experimental data for systems where all-atom simulations are impractical, such as relative folding free energies of protein mutants [5].

The Scientist's Toolkit

Essential Research Reagents and Computational Solutions

Table 3: Key Research Reagents and Computational Tools for Protein Folding Studies

Tool/Solution	Type	Function in Folding Studies	Example Applications
Amber	Software Package	All-atom and coarse-grained MD simulations [22]	Protein folding with ff19SB force field [22]
Gromacs	Software Package	High-performance MD simulations [22]	Martini 3 coarse-grained simulations [22]
CHARMM	Force Field/Software	Molecular mechanics and dynamics [83]	All-atom parameter sets (CHARMM22, CHARMM36) [83] [22]
Martini 3	Coarse-Grained Force Field	Accelerated dynamics of biomolecules [22]	IDP simulations with scaled protein-water interactions [22]
SIRAH	Coarse-Grained Force Field	CG simulations with explicit solvent [22]	IDP and folded protein studies [22]
UNRES	Coarse-Grained Force Field	Physics-based protein folding simulations [22]	Protein folding prediction [22]
CGSchNet	Machine-Learned CG Model	Transferable neural network force field [5]	Extrapolative MD on new sequences [5]
Markov State Models	Analysis Framework	Identifying metastable states and kinetics [84]	Hsp70 chaperone cycling and client folding [84]

Discussion

Interpretation of Comparative Results

The development of machine-learned transferable coarse-grained models represents a significant advancement in protein simulation methodology. These models successfully address the traditional trade-off between computational efficiency and physical accuracy, enabling researchers to study protein folding and dynamics on timescales and system sizes previously inaccessible to all-atom simulations [5].

For small fast-folding proteins, the close agreement between CG and all-atom free energy landscapes demonstrates that the essential physics of folding can be captured at reduced resolution [5]. The ability of CG models to stabilize native-like states with high native contact fractions and low RMSD values indicates that the machine-learned multi-body interactions effectively encode the cooperative effects necessary for specific folding. Furthermore, the reproduction of misfolded states observed in all-atom simulations suggests that these models capture not only the global energy minimum but also the intricate features of the folding landscape [5].

For intrinsically disordered proteins and larger systems, CG models provide particularly valuable insights. The accurate modeling of IDP fluctuations [5] and the ability to simulate folding of proteins beyond the practical limits of all-atom MD [5] significantly expand the scope of addressable research questions in protein biophysics.

Limitations and Considerations

Despite these advances, important limitations remain. The observed challenges with certain structural motifs, such as the mixed α/β structure of BBA where the native state appears as a local rather than global minimum [5], highlight areas for future model refinement. Similarly, the slightly higher sequence fluctuations observed in CG simulations of larger proteins compared to all-atom references [5] suggest potential directions for improving the representation of local flexibility.

When employing coarse-grained models, researchers should consider:

• System Suitability: CG models perform best when the systems of interest share similarities with the training data in terms of sequence features and structural motifs.
• Validation Strategy: Important findings should be validated against available experimental data or targeted all-atom simulations where feasible.
• Force Field Selection: The choice between all-atom and coarse-grained approaches should be guided by research questions, required timescales, and available computational resources [22].

Future Directions

The integration of machine learning with physical modeling continues to advance simulation capabilities. Future developments may include:

• Multi-Scale Frameworks: Seamless integration between all-atom and coarse-grained resolutions for adaptive resolution simulations.
• Enhanced Sampling: Combination of machine-learned CG models with advanced sampling techniques for even greater efficiency.
• Expanded Applicability: Continued broadening of model transferability across diverse protein classes and environmental conditions.

These advancements will further bridge the gap between computational efficiency and physical accuracy, providing researchers with increasingly powerful tools for understanding protein folding and function.

The development of coarse-grained (CG) molecular dynamics parameters for simulating protein folding requires rigorous experimental validation to ensure their predictive accuracy and transferability. This application note details key experimental protocols and benchmarks for validating three critical aspects of protein behavior: folding rates, Φ-values, and mutational effects on stability and binding. The emergence of machine-learned coarse-grained models, which can accelerate dynamics by several orders of magnitude while preserving thermodynamic properties, makes such validation frameworks increasingly essential for researchers and drug development professionals [5] [6]. We frame these methods within the broader context of establishing reliable, computationally efficient alternatives to all-atom molecular dynamics for studying small protein folding.

Core Validated Applications of Coarse-Grained Models

Prediction of Folding Rates and Mechanisms

Protocol 1: Validating Folding Rates and Free Energy Landscapes

• Objective: To assess a CG model's ability to reproduce experimental protein folding rates and correctly identify folding mechanisms (e.g., two-state vs. multi-state) [85].
• Workflow:
  • System Preparation: Initialize CG simulations from both extended and native-state configurations for a set of well-characterized small proteins (e.g., Chignolin, TRP-cage, Fyn SH3 domain).
  • Enhanced Sampling: Perform parallel-tempering or replica-exchange molecular dynamics simulations to ensure converged sampling of the folding/unfolding equilibrium [5].
  • Free Energy Construction: Calculate the free energy landscape as a function of reaction coordinates, typically the fraction of native contacts (Q) and the root-mean-square deviation (RMSD) from the native structure.
  • Rate Calculation: For two-state folders, the folding rate ((kf)) can be estimated from the free energy barrier height ((\Delta G^\ddagger)) using the relation (kf = k0 \exp(-\Delta G^\ddagger / kB T)), where (k0) is a pre-exponential factor.
  • Validation: Compare computed (kf) values and the presence of folding intermediates with experimental data from stopped-flow kinetics or other biophysical measurements [86] [85].

table 1: Comparison of predicted and experimental folding rates for single-domain proteins

Protein (SCOP Class)	Experimental kf (s⁻¹)	CG Model Predicted kf (s⁻¹)	Folding Mechanism
src SH3 (all β)	12 - 18 [85]	~15 [85]	Two-state
En-HD (all α)	~30,000 [85]	~28,000 [85]	Two-state
CI2 (α + β)	~50 [85]	~45 [85]	Two-state
Chignolin	~4,000,000 [6]	Consistent with atomistic reference [5]	Two-state
TRP-cage	~40,000 [6]	Consistent with atomistic reference [5]	Two-state

Φ-Value Analysis for Transition-State Characterization

Protocol 2: Computational Φ-Value Analysis

• Objective: To validate the accuracy of a CG model's predicted folding transition state ensemble (TSE) by comparing computed Φ-values with experimental data [87].
• Background: Φ-value analysis is a key experimental benchmark for characterizing the structure of the TSE. A Φ-value of ~1 indicates the mutated residue is fully structured in the TSE, while a value of ~0 indicates it is unstructured.
• Workflow:
  • TSE Identification: From free energy landscapes (see Protocol 1), identify configurations belonging to the transition state ensemble (typically defined as structures with Q around 0.5-0.6).
  • In Silico Mutagenesis: For a set of residues with experimental Φ-values, computationally introduce the same point mutations into all members of the TSE and the native state (N).
  • Free Energy Change Calculation: Calculate the change in free energy of the TSE ((\Delta \Delta G^\ddagger)) and the native state ((\Delta \Delta GN)) due to the mutation. The Φ-value is computed as: [ \Phi = \frac{\Delta \Delta G^\ddagger}{\Delta \Delta GN} ]
  • Precision Considerations: Note that the precision of Φ-values is highly sensitive to experimental error when (\Delta \Delta G_N) is small (typically < 7 kJ/mol) [86]. Validation requires comparison across multiple mutations with significant stability changes.

table 2: Key considerations for experimental and computational φ-value determination

Aspect	Considerations & Challenges	Recommendations
Precision	Poor when ΔΔG < 7 kJ/mol; subject to inter-laboratory variability [86]	Use rates at non-zero denaturant or assume invariant chevron slopes (mf, mu) to improve precision [86].
Data Quality	Scatter in φ-values can stem from experimental error or genuine mechanistic changes [86].	Use blind, replicate measurements to assess true precision.
Computational Correlation	Calculated ΔΔG of the TSE can be correlated with experimental φ-values [87].	Use the method to validate all-atom or CG models of the TSE.

Predicting Mutational Effects on Stability and Binding

Protocol 3: Predicting Binding Energy Changes via Folding Stability

• Objective: To predict the effect of mutations on protein-protein binding affinity ((\Delta \Delta G_{bind})) using a transfer-learning approach based on folding energies [88] [89].
• Theoretical Basis: The binding free energy can be parameterized as a difference in folding stabilities, leveraging the state function property of free energy [88]: [ \Delta \Delta G{bind}(A:B \rightarrow A':B) = \Delta \Delta G{fold}(A:B \rightarrow A':B) - \Delta \Delta G_{fold}(A \rightarrow A') ]
• Workflow (StaB-ddG Method):
  • Initialization: Use a pre-trained inverse folding model (e.g., ProteinMPNN) as a proxy for sequence log-likelihood, which correlates with folding stability.
  • Fine-Tuning: Fine-tune the model in two stages:
    • First, on high-throughput experimental data of folding stability changes ((\Delta \Delta G{fold})).
    • Second, on more limited datasets of binding energy changes ((\Delta \Delta G{bind})).
  • Prediction & Validation: For a mutant complex, predict the (\Delta \Delta G_{bind}) using the fine-tuned model (StaB-ddG) and validate against experimental measurements or expensive force-field calculations (e.g., FoldX). This approach can match the accuracy of FoldX while being over 1,000 times faster [88].

table 3: Performance of methods for predicting mutational effects on binding energy

Method	Principle	Speed	Accuracy (vs. Experiment)
StaB-ddG	Deep learning; transfer learning from folding stability [88] [89]	Very Fast (~1000x Faster than FoldX)	Matches state-of-the-art (FoldX)
FoldX (Empirical Force Field)	Empirical force field based on full atomic details [88]	Slow	State-of-the-art
Other Deep Learning Predictors	Direct learning from limited binding data [88]	Fast	Underperforms due to data scarcity

Experimental Workflow and Logical Framework

The following diagram illustrates the integrated workflow for the experimental validation of coarse-grained models, connecting the protocols described above.

Workflow for CG Model Validation

The Scientist's Toolkit: Essential Research Reagents and Materials

table 4: Key research reagents and computational tools for validation experiments

Item / Resource	Function / Description	Application in Validation
Stopped-Flow Spectrometer	Instrument for rapid mixing and monitoring of fast folding/unfolding kinetics (millisecond resolution) [86].	Protocol 1: Measuring experimental folding rates (kf) for benchmark proteins.
Denaturants (e.g., Urea, GdmCl)	Chemical denaturants used to unfold proteins and map chevron plots (ln(k) vs. [denaturant]) [86].	Protocol 1 & 2: Determining folding/unfolding rates and m-values; essential for φ-value analysis.
Curated Protein Folding/Binding Datasets	Publicly available datasets of experimental folding (ΔΔGfold) and binding (ΔΔGbind) energies for mutants [88].	Protocol 3: Training and benchmarking computational predictors like StaB-ddG.
Coarse-Grained Simulation Software (e.g., CGSchNet)	Machine-learning based MD software implementing neural network CG force fields [5] [6] [27].	All Protocols: Running the CG simulations that generate predictions for validation.
ProteinMPNN	A pre-trained inverse folding neural network model for protein sequence design [88].	Protocol 3: Serves as a foundation model for initializing zero-shot predictors of folding stability.
FoldX	An empirical force field method for predicting the effect of mutations on stability and interactions [88].	Protocol 3: Used as a state-of-the-art benchmark for comparing the accuracy of new predictors.

The experimental protocols and benchmarks outlined herein provide a robust framework for validating coarse-grained molecular dynamics parameters. The ability of modern, machine-learned CG models to predict quantitative experimental observables—such as folding rates, Φ-values, and mutational effects—with high accuracy while offering computational speedups of several orders of magnitude, demonstrates their growing maturity and value [5] [6]. This validation is crucial for establishing these models as reliable tools for fundamental research in protein folding and for applied drug development efforts, such as predicting the functional consequences of mutations in therapeutic proteins.

Molecular dynamics (MD) simulations serve as a cornerstone of atomistic-scale analysis across numerous disciplines, including drug discovery and structural biology [90] [91]. For studies focused on small protein folding, a fundamental challenge has been the prohibitive computational cost of achieving ab initio accuracy, which is essential for reliably capturing phenomena like folding pathways and native state stability [90] [92]. Current state-of-the-art simulations using ab initio MD (AIMD) or machine-learning MD (MLMD) on general-purpose central/graphics processing units (CPU/GPU) are intrinsically limited by "memory wall" and "power wall" bottlenecks [90]. This imposes severe restrictions on the feasible simulation size and timescales for protein folding research.

Within this context, coarse-grained (CG) models present a strategic alternative for studying protein folding by reducing the number of explicit degrees of freedom [16] [7]. The physical basis of coarse-grained molecular dynamics stems from the potential of mean force, where the motion of CG sites is governed by this potential along with friction and stochastic forces resulting from integrating out the secondary, fine-grained degrees of freedom [16]. This review analyzes a paradigm-shifting approach—the special-purpose MD processing unit (MDPU)—which promises to mitigate traditional hardware limitations. We provide a quantitative analysis demonstrating computational cost reductions of 3-4 orders of magnitude, detail protocols for leveraging these advances in protein folding studies, and contextualize these developments within the broader framework of coarse-grained MD parameters.

Results & Performance Analysis

Quantitative Analysis of Speed and Power Efficiency

The performance gap between traditional computational platforms and emerging specialized hardware is profound. The table below summarizes key performance metrics for different MD approaches, highlighting the dramatic efficiency gains offered by the MDPU architecture.

Table 1: Performance Metrics of MD Simulation Approaches

Processing Platform	Time Consumption (ηt, sec/step/atom)	Power Consumption (ηp, Watt)	Accuracy (Energy RMSE, meV/atom)	Key Methodological Basis
Ab Initio MD (AIMD) on CPU/GPU	~10⁻⁶ [90]	Up to 10¹ MW (large clusters) [90]	N/A (Reference Accuracy)	Density Functional Theory [90]
Machine-Learning MD (MLMD) on CPU/GPU	~10⁻⁹ [90]	High (requires extensive parallelization) [90]	~1-10 meV/atom [90]	Machine-learned interatomic potentials [90]
Special-Purpose MDPU	~10⁻¹² (109x vs AIMD) [90]	Dramatically reduced (103x vs MLMD) [90]	~1.66-85.35 meV/atom [90]	Computing-in-Memory (CIM) architecture; Co-design of algorithm/hardware/software [90]
Coarse-Grained MD (CGMD)	Varies by resolution; significantly faster than all-atom [16]	Lower due to reduced particle count	Varies by force field parameterization [16] [7]	Potential of Mean Force; Langevin dynamics [16]

The MDPU achieves these performance breakthroughs through a fundamental rethinking of computational architecture. By implementing a computing-in-memory (CIM) engine, it directly mitigates the von Neumann "memory wall" bottleneck, where traditional CPUs/GPUs spend over 90% of their time and power shuttling data rather than performing calculations [90]. This hardware innovation is complemented by algorithmic replacements of heavy-duty calculations with lightweight equivalent operations derived from revised machine-learning approaches [90].

Accuracy Validation Across Biomolecular Systems

Despite the dramatic speed improvements, the MDPU approach maintains quantum mechanical accuracy, which is critical for reliable protein folding studies. Validation across diverse systems demonstrates its capability to capture subtle energetic phenomena:

Table 2: Accuracy Benchmarks Across Various Molecular Systems

Molecular System	Energy Prediction Error (εe, meV/atom)	Force Prediction Error (εf, meV/Å)	Relevance to Protein Folding Studies
Au (Elemental)	85.35	173.20	Validation under extreme conditions [90]
Cu (Elemental)	1.84	16.44	Mechanical property prediction [90]
Mg (Elemental)	2.47	13.91	Stacking fault energy calculation [90]
H₂O (Binary)	7.62	148.24	Hydrogen bonding & solvation effects [90]
GeTe (Binary)	2.95	110.69	Phase change & structural transitions [90]
Li₁₀GeP₂S₁₂ (Quaternary)	1.66	59.67	Complex biomolecular analog [90]

The accuracy demonstrated across these systems, particularly the ability to capture hydrogen bonding in water (εe = 7.62 meV/atom) and subtle structural transitions in complex systems like Li₁₀GeP₂S₁₂ (εe = 1.66 meV/atom), establishes the MDPU's capability for studying protein folding energetics [90]. These errors are significantly smaller than the typical interatomic bond energy (10⁰-10¹ eV/atom) and approach the default convergence thresholds for force in popular DFT packages (23.1-40.0 meV/Å) [90].

Experimental & Computational Protocols

Workflow for High-Speed Protein Folding Studies

The following diagram illustrates the integrated workflow combining coarse-grained modeling principles with specialized hardware acceleration for protein folding studies:

Protocol 1: Coarse-Grained Structure Modeling for MDPU Simulation

Objective: Generate accurate initial structures for MDPU-accelerated folding studies using coarse-grained representations.

• Amino Acid Sequence Preparation

  • Obtain protein sequence in FASTA format from databases (e.g., UniProt, NCBI) [93].
  • Perform multiple sequence alignment to identify evolutionary constraints [92].

• Secondary Structure Prediction

  • Utilize PSIPRED or similar tools to predict secondary structure elements (α-helices, β-sheets, coils) [93].
  • Use predictions to inform coarse-grained mapping strategy.

• Coarse-Grained Mapping Selection

  • For backbone-dominated studies: Implement Cα-based mapping (1 site per residue) [7].
  • For larger systems: Apply essential dynamics coarse-graining (ED-CG) to identify dynamic domains [7].
  • Define resolution level based on research question (e.g., 4:1 to 10:1 heavy atom to CG site ratio) [16].

• Tertiary Structure Modeling (Choose one approach)

  • De Novo Modeling (Template-Free):
    • AlphaFold2: Use Colab notebook for prediction without templates [93] [92].
    • Robetta: Submit sequence to RoseTTAFold server for model generation [93].
    • trRosetta: Utilize web platform for distance-based structure prediction [93].
  • Template-Based Modeling:
    • MOE Homology Modeling: Perform BLAST search to identify templates with >40% sequence identity [93].
    • I-TASSER: Use automated platform for template identification and threading [93].

• Model Quality Assessment

  • Calculate phi-psi plots to evaluate backbone torsion angles [93].
  • Use ERRAT for overall quality assessment [93].
  • Compare B-factors with experimental data if available [7].

Protocol 2: MDPU-Accelerated Folding Dynamics

Objective: Simulate protein folding pathways with ab initio accuracy at significantly reduced computational cost.

• System Preparation

  • Convert CG structure to MDPU-compatible format.
  • Apply solvation model appropriate for CG representation [16].
  • Set ionic concentration to physiological levels (e.g., 150mM NaCl).

• Force Field Parameterization

  • Implement potential of mean force derived from all-atom simulations [16].
  • Parameterize friction coefficients (Γ) for Langevin dynamics [16].
  • Set stochastic force terms to maintain constant temperature [16].

• MDPU Simulation Configuration

  • Initialize with computing-in-memory architecture to minimize data transfer [90].
  • Configure pipeline processing optimized for FPGA/ASIC hardware [90].
  • Set neural network inference for energy and force calculations [90].

• Production Simulation

  • Equilibrate system with position restraints on protein backbone.
  • Gradually remove restraints over 1-5ns simulation time.
  • Run production simulation for sufficient duration to observe folding/unfolding events.
  • Employ enhanced sampling techniques if needed for rare events.

• Trajectory Analysis

  • Monitor root mean square deviation (RMSD) of backbone atoms [93].
  • Calculate radius of gyration (Rg) to track compaction [93].
  • Identify native contacts formation using contact maps.
  • Analyze free energy landscape as function of collective variables.

The Scientist's Toolkit: Research Reagent Solutions

Table 3: Essential Computational Tools for CG-MD Protein Folding Studies

Tool/Category	Specific Examples	Function & Application	Hardware Compatibility
Structure Prediction	AlphaFold2, Robetta, trRosetta [93]	De novo protein structure prediction from sequence	CPU/GPU, Cloud-based [92]
Homology Modeling	MOE, I-TASSER [93]	Template-based structure construction	CPU/GPU [93]
Coarse-Grained Force Fields	MARTINI, UNICORN [16]	Define interactions between CG sites using potential of mean force	CPU/GPU, MDPU [16]
Specialized MD Hardware	MDPU, ANTON [90]	Hardware-accelerated molecular dynamics	Dedicated systems [90]
Dynamic Network Models	DENM, CG-DENM [7]	Elastic network models for large-scale conformational changes	CPU/GPU [7]
Analysis & Visualization	MOE, VMD, PyMOL [93]	Trajectory analysis, visualization, and property calculation	CPU/GPU [94]

Discussion & Future Perspectives

The integration of specialized hardware like MDPU with advanced coarse-grained methodologies represents a transformative development for protein folding studies. The demonstrated computational cost reductions of 3-4 orders of magnitude directly address the critical barrier of simulation timescales, potentially enabling researchers to access biologically relevant folding times that were previously computationally prohibitive [90].

The physical basis of coarse-grained MD—where dynamics emerge from the potential of mean force with friction and stochastic forces accounting for integrated degrees of freedom—provides a rigorous foundation for these accelerated simulations [16]. When combined with the MDPU architecture that mitigates memory wall bottlenecks through computing-in-memory technology, this approach offers a pathway to maintain quantum mechanical accuracy while achieving unprecedented simulation speeds [90].

For the specific application of small protein folding studies, these advances enable several new research directions:

• High-throughput folding landscape mapping for mutational studies and protein design applications.
• Direct simulation of folding mechanisms for proteins that lack experimental structural data.
• Multi-scale approaches that combine CG folding simulations with all-atom refinement of populated states.

Future development efforts should focus on optimizing CG force fields specifically for MDPU architectures, improving the integration between deep learning structure prediction and dynamics simulation, and validating folding pathways against emerging experimental techniques. As these methods mature, the combination of algorithmic innovations like differentiable deep learning [92] with specialized hardware promises to further bridge the gap between computational feasibility and biological reality in protein folding research.

Coarse-grained (CG) molecular dynamics has become an indispensable tool for simulating biological processes, such as protein folding, that occur at time and length scales beyond the reach of all-atom (AA) simulations [16] [3]. By representing groups of atoms as single interaction sites, or "beads," CG models dramatically reduce the number of degrees of freedom, enabling simulations that are orders of magnitude faster and covering vastly larger scales [9]. Despite these advantages, CG models face profound challenges that limit their predictive accuracy and transferability. This application note details the central limitations of current CG models in the context of small protein folding studies, providing a critical analysis for researchers and drug development professionals. It further outlines established protocols for model evaluation and discusses emerging frontiers, particularly machine learning (ML), that aim to overcome these hurdles.

Core Limitations of Contemporary CG Models

The development of a truly universal and predictive CG model for proteins remains a formidable challenge. The primary obstacles can be categorized into several key areas, summarized in Table 1 below.

Table 1: Key Limitations in Coarse-Grained Protein Folding Models

Limitation Category	Specific Challenge	Impact on Protein Folding Studies
Transferability & Chemical Specificity	Limited accuracy on sequences not used for parameterization; difficulty capturing effects of point mutations [5] [9].	Hinders predictive simulation of novel protein sequences or mutants; reduces utility in drug discovery.
Multi-Body Interactions & Implicit Solvation	Inefficient modeling of essential multi-body interaction terms for correct protein thermodynamics [5].	Can lead to incorrect stabilization of non-native states or inaccurate folding free energy landscapes [1].
Dynamics & Timescales	Smoothed energy landscapes accelerate dynamics but break direct link to real-time kinetics [1] [9].	Prevents accurate prediction of folding rates and pathways; relationship to AA/real time is often ambiguous [16].
Force Field Parameterization	Dependence on AA simulations for parameterization propagates inaccuracies from the reference force field [6].	Model is only as good as its reference data; errors in the AA force field become embedded in the CG model.
Spatial Resolution & Structural Accuracy	Loss of atomic detail can miss key interactions; instability in folded states of larger proteins [5].	Limits application to systems where specific side-chain interactions or precise geometry are critical.

The Challenge of Transferability and Chemical Specificity

A fundamental limitation of many CG models is their lack of transferability—the ability to perform accurately on protein sequences not included in their training set [9]. For instance, while some models succeed on proteins with high sequence similarity to training data, they struggle with sequences sharing low (<40%) similarity [5]. This directly impacts the study of protein mutants, a crucial area in understanding disease and drug design. A related issue is the loss of chemical specificity, where simplified representations fail to adequately capture the unique physicochemical properties of different amino acid side chains, affecting the model's ability to distinguish between closely related proteins or the effects of single-point mutations [9].

The Many-Body Interaction Problem

From a statistical mechanics perspective, the correct potential for a CG model is the Potential of Mean Force (PMF), which inherently includes complex multi-body correlations [16]. However, many popular CG force fields rely on simple, pairwise interaction potentials copied from all-atom force fields. This simplification is a major hindrance to accurately representing protein thermodynamics, as it cannot fully capture the intricate interactions that stabilize native protein structures, particularly the effects of implicit solvation [5]. For example, even advanced models like UNRES and AWSEM have been noted to sometimes fail to capture alternative metastable states or accurately model the stability of proteins with mixed α/β motifs [5].

Dynamics, Kinetics, and the Timescale Dilemma

While CG simulations famously accelerate dynamics by smoothing the energy landscape, this creates a problem for kinetic studies. The motion of CG sites is governed by the PMF and friction/stochastic forces from the integrated atomic degrees of freedom, making Langevin dynamics a natural choice [16]. However, the simplified landscape distorts the true energy barriers, breaking a direct correspondence between simulation time and real time. Although this allows faster observation of folding events, it makes the prediction of experimentally measurable folding rates unreliable [1]. Furthermore, the common representation of CG sites as spherical point masses ignores molecular anisotropy, which can adversely affect the accuracy of simulated dynamics [16].

Experimental Protocols for Evaluating CG Model Performance

To systematically assess the limitations of a CG model, a standardized set of validation protocols is essential. The following workflow, detailed in Table 2, outlines key experiments cited in recent literature.

Diagram 1: Workflow for CG model evaluation. This workflow outlines the key steps for a comprehensive assessment of a coarse-grained model's performance.

Table 2: Protocol for Benchmarking CG Models in Protein Folding

Protocol Step	Description	Key Metrics & Outputs
1. Free Energy Landscape Calculation	Use enhanced sampling methods (e.g., Parallel Tempering, Metadynamics) to construct a converged free energy surface as a function of reaction coordinates like fraction of native contacts (Q) and Cα root-mean-square deviation (RMSD) [5] [6].	- Number and stability of metastable states (folded, unfolded, intermediates).- Relative free energy differences between states.- Structural characterization of minima.
2. Folding/Unfolding Simulations	Initiate long, constant-temperature simulations from both extended and native states. Compare folding mechanisms and transition paths to all-atom reference simulations where available [5] [1].	- Observation of multiple folding/unfolding events.- Folding pathways and mechanisms (e.g., order of structure formation).- Identification of potential misfolded traps.
3. Transferability Testing	Simulate proteins with low sequence similarity (<40%) to any protein in the model's training set. A critical test is the prediction of relative folding free energies for point mutants [5] [95].	- Accuracy in predicting native structure for novel sequences.- Correlation of predicted ΔΔGfold with experimental data.- Structural fluctuations (Cα RMSF) compared to AA or experiment.
4. Dynamic Property Validation	Analyze internal dynamics and fluctuations within the folded state. Compare to experimental data (e.g., NMR relaxation) or all-atom simulations [5] [96].	- Cα root-mean-square fluctuations (RMSF) per residue.- Comparison of domain motions and flexibility patterns.

Detailed Protocol: Free Energy Landscape Analysis

Objective: To determine the thermodynamic ensemble and identify all relevant metastable states in the protein's conformational landscape.

Methodology:

• System Setup: Place the CG protein in a suitable simulation box. For implicit solvent models, no explicit solvent is needed. For explicit solvent CG models (e.g., MARTINI), use the appropriate CG water and ion parameters.
• Reaction Coordinates: Define collective variables (CVs). The most common and effective CVs for folding are:
  • Fraction of Native Contacts (Q): Calculated as the number of CG native contacts present in a given configuration divided by the total number in the native state. A contact can be defined based on a distance cutoff in the reference structure.
  • Cα Root-Mean-Square Deviation (RMSD): Measures the global structural deviation from the native state.
• Enhanced Sampling: Employ Parallel Tempering (Replica Exchange) molecular dynamics. This involves running multiple simulations (replicas) of the same system at different temperatures (or with different Hamiltonian parameters). Exchanges between replicas are attempted periodically based on the Metropolis criterion, which ensures efficient sampling over energy barriers.
• Analysis:
  • Use the weighted histogram analysis method (WHAM) to combine data from all replicas and calculate the free energy as a function of the chosen CVs: ( G(Q, RMSD) = -k_B T \, \text{ln} \, P(Q, RMSD) ), where ( P ) is the probability distribution.
  • Identify free energy minima (stable states) and the saddle points (energy barriers) between them.

The Scientist's Toolkit: Research Reagent Solutions

Table 3: Essential Resources for CG Protein Folding Research

Resource	Function in Research	Example Use Case
CGSchNet Model [5]	A machine-learned, transferable CG force field for proteins. Uses a neural network to model many-body interactions from AA data.	Extrapolative MD on new protein sequences not used in training; predicting folding landscapes and mutant effects.
Gō-MARTINI Model [96]	A hybrid model combining the MARTINI force field with a structure-based Gō model. Good for protein-membrane systems.	Simulating large-scale structural dynamics of proteins on lipid membranes, such as in membrane remodeling.
UNRES Force Field [16] [5]	A physics-based CG force field developed for protein folding studies.	Studying protein folding pathways and the thermodynamics of small, fast-folding proteins.
HMDM Benchmark Dataset [95]	A benchmark dataset of high-quality homology models for evaluating Model Quality Assessment (MQA) methods.	Testing the performance of CG models or scoring functions on accurate, single- and multi-domain protein structures.
Multi-millisecond AA MD Dataset [6]	A unique dataset of ~9 ms of unbiased all-atom MD simulations for 12 different proteins.	Serving as a reference for training and validating bottom-up CG models via force-matching or relative entropy minimization.

Frontiers and Future Directions

The most promising frontier in overcoming current limitations is the integration of machine learning (ML). ML approaches, particularly deep learning, are being used to create CG potentials that can faithfully represent the many-body PMF from large-scale AA simulation data [14] [5] [6]. For example, the CGSchNet model demonstrates that a single, sequence-transferable potential can predict folded states, fluctuations of disordered proteins, and relative folding free energies for mutants, rivaling the qualitative accuracy of AA-MD at a fraction of the computational cost [5]. These methods leverage techniques like variational force-matching to learn effective forces directly from atomistic data, promising a new generation of CG models that are both highly accurate and broadly transferable [6]. The continued development and benchmarking of these models, using rigorous protocols as outlined above, represent the critical path forward for the field of coarse-grained protein simulation.

Conclusion

The integration of machine learning with coarse-grained molecular dynamics represents a paradigm shift in protein folding studies, enabling accurate prediction of folding landscapes, mechanisms, and thermodynamics at a fraction of the computational cost of all-atom simulations. The development of transferable, chemically aware CG force fields through deep learning and advanced optimization techniques has overcome traditional limitations of system-specific parameterization. While challenges remain in capturing complex multi-body interactions and ensuring robust physical understanding, these advances open new frontiers for simulating larger biomolecular systems, protein-ligand interactions, and pathogenic mutations relevant to disease. As CG methodologies continue evolving toward greater accuracy and physical realism, they promise to accelerate drug discovery by providing unprecedented insights into protein folding diseases and enabling high-throughput virtual screening of therapeutic interventions at biologically relevant timescales.

References

• 1. Challenges in protein folding simulations: Timescale, ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC3032381/]
• 2. Successes and challenges in simulating the folding of large ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC6952611/]
• 3. Recent Advances in Coarse-Grained Models for ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC6696403/]
• 4. Insights from Coarse-Grained Gō Models for Protein Folding ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC2672008/]
• 5. Navigating protein landscapes with a machine-learned ... [https://www.nature.com/articles/s41557-025-01874-0]
• 6. Machine learning coarse-grained potentials of protein ... [https://www.nature.com/articles/s41467-023-41343-1]
• 7. COARSE-GRAINED MODELING OF PROTEIN ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC4230303/]
• 8. Understanding Protein Solution Phase Behavior via ... [https://www.nist.gov/programs-projects/understanding-protein-solution-phase-behavior-coarse-grained-modeling]
• 9. Recent Advances in Transferable Coarse-Grained ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC5366245/]
• 10. OpenCafeMol: A coarse-grained biomolecular simulator on ... [https://www.sciencedirect.com/science/article/pii/S0006349525004473]
• 11. Zwitterionic lipid assemblies: Molecular dynamics studies of ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC2876730/]
• 12. Protein–ligand binding with the coarse-grained Martini model [https://www.nature.com/articles/s41467-020-17437-5]
• 13. A Standardized Benchmark for Machine-Learned ... [https://arxiv.org/html/2510.17187v1]
• 14. Recent advances in machine learning and coarse-grained ... [https://www.sciencedirect.com/science/article/pii/S0006349525003790]
• 15. : An Overview | MRS Online Proceedings Library Energy Landscapes [https://link.springer.com/article/10.1557/PROC-700-S8.5]
• 16. Theory and Practice of Coarse-Grained Molecular ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC8465211/]
• 17. Genomic Energy - PMC Landscapes [https://pmc.ncbi.nlm.nih.gov/articles/PMC5300775/]
• 18. Coarse-Grained Simulations of Protein-Protein Association [https://pmc.ncbi.nlm.nih.gov/articles/PMC3443792/]
• 19. Characterizing Protein Energy Landscape by Self-Learning ... [https://www.sciencedirect.com/science/article/pii/S0006349510010362]
• 20. Protocol for the development of coarse-grained structures ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC10399882/]
• 21. GENESIS CGDYN: large-scale coarse-grained MD ... [https://www.nature.com/articles/s41467-024-47654-1]
• 22. A Practical Guide to All-Atom and Coarse-Grained ... [https://www.mdpi.com/1422-0067/25/12/6698]
• 23. Implicit-solvent coarse-grained modeling for polymer ... [https://pubs.rsc.org/en/content/articlelanding/2019/sm/c9sm01211g]
• 24. Bottom-Up Coarse-Grained Modeling of DNA [https://www.frontiersin.org/journals/molecular-biosciences/articles/10.3389/fmolb.2021.645527/full]
• 25. A Review of Physics-Based Coarse-Grained Potentials for ... [https://www.sciencedirect.com/science/article/abs/pii/B9780444594402000053]
• 26. Investigation of protein folding by coarse-grained molecular ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC2849147/]
• 27. Machine Learning of Coarse-Grained Molecular Dynamics ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC6535777/]
• 28. A beginner's approach to deep learning applied to VS and MD ... [https://jcheminf.biomedcentral.com/articles/10.1186/s13321-025-00985-7]
• 29. Coarse graining molecular dynamics with graph neural networks [https://pmc.ncbi.nlm.nih.gov/articles/PMC7671749/]
• 30. Coarse-Graining with Equivariant Neural : A Path Toward... Networks [https://pmc.ncbi.nlm.nih.gov/articles/PMC10726966/]
• 31. NeurIPS Poster Efficiently incorporating quintuple interactions into... [https://neurips.cc/virtual/2023/poster/71146]
• 32. Peering inside the black box by learning the relevance of many - body ... [https://www.nature.com/articles/s41467-025-65863-0?error=cookies_not_supported&code=57bc09e4-c69f-47ed-a42f-8f96b5159769]
• 33. Folding of Small Proteins Using Constrained Molecular ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC3124276/]
• 34. Introducing a force-matched united atom force field to ... [https://www.sciencedirect.com/science/article/pii/S026412752400203X]
• 35. Deep generative modeling of temperature-dependent ... - PMC [https://pmc.ncbi.nlm.nih.gov/articles/PMC12627665/]
• 36. Reweighting configurations generated by transferable ... [https://pubs.rsc.org/en/content/articlehtml/2025/me/d4me00198b]
• 37. Refining coarse-grained molecular topologies: a Bayesian ... [https://www.nature.com/articles/s41524-025-01729-9]
• 38. Parametrization of a new small molecule [https://cgmartini.nl/docs/tutorials/Martini3/Small_Molecule_Parametrization/]
• 39. [2501.02707] Refining Coarse-Grained Molecular Topologies [https://arxiv.org/abs/2501.02707]
• 40. A primer on Bayesian optimisation - NUBO [https://nubopy.com/_build/html/bayesian_optimisation.html]
• 41. Scalable Bayesian Optimization for High-Dimensional ... [https://arxiv.org/html/2506.22533v1]
• 42. Bayesian optimization in complex scenarios: possible! (BO ... [https://medium.com/mlearning-ai/bayesian-optimization-in-complex-scenarios-possible-bo-series-ii-443729730c13]
• 43. Understanding protein folding: Small proteins in silico - PMC [https://pmc.ncbi.nlm.nih.gov/articles/PMC2244683/]
• 44. SimpleFold: Folding Proteins is Simpler than You Think [https://arxiv.org/html/2509.18480v2]
• 45. Sampling of structure and sequence space of small protein ... [https://www.nature.com/articles/s41467-022-34937-8]
• 46. Folding Trp-Cage to NMR Resolution Native Structure Using a ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC1304993/]
• 47. Folding Trp-cage to NMR resolution native structure using ... [https://arxiv.org/abs/physics/0405134]
• 48. Accelerated Molecular Dynamics Simulations of Protein Folding [https://pmc.ncbi.nlm.nih.gov/articles/PMC4487363/]
• 49. Differentiable molecular simulation can learn all the ... [https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0256990]
• 50. The Unfolded State of the Villin Headpiece Helical Subdomain [https://pmc.ncbi.nlm.nih.gov/articles/PMC4805113/]
• 51. Monte Carlo, harmonic approximation, and coarse - graining ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC2924683/]
• 52. Enhanced Sampling Techniques in Molecular Dynamics ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC4339458/]
• 53. A Deep Graph Network–Enhanced Sampling Approach to ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC8116519/]
• 54. Enhanced Sampling for Efficient Learning of Coarse ... [https://arxiv.org/abs/2510.11148]
• 55. Investigating whether deep learning models for co-folding ... [https://www.nature.com/articles/s41467-025-63947-5]
• 56. Force Field Bias in Protein Folding Simulations - PMC [https://pmc.ncbi.nlm.nih.gov/articles/PMC2711430/]
• 57. Grappa – a machine learned molecular mechanics force field [https://pmc.ncbi.nlm.nih.gov/articles/PMC11734696/]
• 58. Quantifying the accuracy of protein simulation models [https://foldingathome.org/2025/09/04/quantifying-the-accuracy-of-protein-simulation-models/]
• 59. How to prevent overfitting — Regularization of ML models [https://medium.com/@ckekula/regularization-of-machine-learning-models-a-mathematical-guide-part-1-4a6bc0b4eef2]
• 60. Overfitting and Regularization in ML [https://www.geeksforgeeks.org/machine-learning/overfitting-and-regularization-in-ml/]
• 61. The Multiple Time-Stepping Method for 3-Body Interactions ... [https://arxiv.org/html/2507.11172v1]
• 62. Peering inside the black box by learning the relevance of ... [https://www.nature.com/articles/s41467-025-65863-0]
• 63. Using Short Molecular Dynamics Simulations to Determine ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC11878205/]
• 64. Intrinsic α-helical and β-sheet conformational preferences [https://pmc.ncbi.nlm.nih.gov/articles/PMC4088981/]
• 65. New class of protein misfolding simulated in high definition [https://science.psu.edu/news/Obrien8-2025]
• 66. Ab initio characterization of protein molecular dynamics ... [https://www.nature.com/articles/s41586-024-08127-z]
• 67. Learning meaningful representations of protein sequences [https://www.nature.com/articles/s41467-022-29443-w]
• 68. Advances in coarse-grained modeling of macromolecular ... [https://www.sciencedirect.com/science/article/abs/pii/S0959440X18300939]
• 69. Operator Forces For Coarse-Grained Molecular Dynamics [https://arxiv.org/html/2506.19628v1]
• 70. A coarse‐grained model for disordered proteins under ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12268113/]
• 71. towards solving the protein-folding problem computationally [https://pmc.ncbi.nlm.nih.gov/articles/PMC6790072/]
• 72. Improving the Physical Realism and Structural Accuracy of ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC3218324/]
• 73. How Robust Are Protein Folding Simulations with Respect ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC3149239/]
• 74. Optimized protein-water interactions and torsional ... [https://www.nature.com/articles/s41467-025-65603-4]
• 75. Comparison of Sequence-Based and Structure ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC1305466/]
• 76. High Energy Channeling and Malleable Transition States [https://www.mdpi.com/2218-273X/12/7/940]
• 77. Uncertainty driven active learning of coarse grained free ... [https://www.nature.com/articles/s41524-023-01183-5]
• 78. Obtaining protein foldability information from computational ... [https://www.sciencedirect.com/science/article/pii/S2001037022003683]
• 79. Native contacts determine protein folding mechanisms in ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC3816414/]
• 80. Navigating protein landscapes with a machine-learned ... [https://communities.springernature.com/posts/navigating-protein-landscapes-with-a-machine-learned-transferable-coarse-grained-model]
• 81. Generation of protein dynamics by machine learning [https://www.sciencedirect.com/science/article/abs/pii/S0959440X25001332]
• 82. Grappa – a machine learned molecular mechanics force field [https://pubs.rsc.org/en/content/articlehtml/2025/sc/d4sc05465b]
• 83. All-Atom and Coarse-Grained Molecular Dynamics ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC3214636/]
• 84. Role of Hsp70 Chaperone in Client-Protein Folding ... [https://pubmed.ncbi.nlm.nih.gov/40873036/]
• 85. Accurate prediction of protein folding mechanisms by simple ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC10587348/]
• 86. On the precision of experimentally determined protein ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC2249776/]
• 87. Calculation of Mutational Free Energy Changes in ... [https://www.sciencedirect.com/science/article/pii/S0006349503745561]
• 88. Predicting mutational effects on protein binding from ... [https://arxiv.org/html/2507.05502v1]
• 89. Predicting mutational effects on protein binding from ... [https://openreview.net/forum?id=En0vq4ICkW¬eId=nrqGPI5pmv]
• 90. High-speed and low-power molecular dynamics ... [https://www.nature.com/articles/s41524-024-01422-3]
• 91. Molecular Dynamics Simulation for All [https://www.sciencedirect.com/science/article/pii/S0896627318306846]
• 92. Folding Revolution [https://hms.harvard.edu/news/folding-revolution]
• 93. Comparison, Analysis, and Molecular Dynamics ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC10525864/]
• 94. From complex data to clear insights: visualizing molecular ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC11043564/]
• 95. A Benchmark Dataset for Evaluating Practical Performance ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC8945737/]
• 96. Optimizing Gō-MARTINI Coarse-Grained Model for F-BAR ... [https://www.frontiersin.org/journals/molecular-biosciences/articles/10.3389/fmolb.2021.619381/full]