Protein Structure and Function: Applications of Bioinformatics Methods - John Rigden 2014

Protein Dynamics: From Structure to Function
Molecular Dynamics Calculations
Limitations and Enhanced Sampling Algorithms

Although Structure/8.html">Molecular Dynamics (MD) simulations have become an integral part of structural biology and have repeatedly provided invaluable insights into biological processes at the atomic level, the field still faces both methodological and computational limitations. Methodological constraints stem from the classical description of atoms and the approximation of atomic interactions using simple potential energy terms instead of the Schrödinger equation. Consequently, Chemical Reactions—such as the breaking and formation of chemical bonds—cannot be described. Furthermore, polarization effects and proton tunneling lie beyond The Scope of classical MD simulations.

The second category of limitations arises from the computational demands of MD simulations. Although bonds are typically treated as spatial constraints imposed on atoms to eliminate the highest-frequency motions, the time step in MD simulations generally cannot exceed 2 fs. A 1-nanosecond simulation, therefore, requires 500,000 force evaluations and an equal number of integration steps. With current algorithms and computing power, timescales on the order of 100 ns are achievable after 3–4 weeks of calculation for a solvated protein with 200 residues.

However, biologically crucial protein motions, such as large-scale conformational transitions, folding, or Denaturation, occur on timescales ranging from microseconds to milliseconds. It thus becomes evident that despite the ever-increasing performance of computing hardware—which has grown roughly a hundredfold every decade—MD simulations will not be able to solve the sampling problem in the foreseeable future solely through increases in computer speed. Therefore, alternative Methods, partly based on MD, have been proposed specifically to address The problem of conformational sampling and the prediction of functionally significant protein motions.

One approach involves reducing the number of particles. Since Proteins are typically studied in a solvent environment, and the majority of the modeled system consists of Water molecules, developing implicit solvent models is a promising way to reduce computational resource requirements (Still et al. 1990; Gosh et al. 1998; Jean-Charles et al. 1991; Luo et al. 2002). Another way to decrease the number of particles is to use so-called coarse-grained models (Bond et al. 2007), in which atoms are grouped into clusters called pseudoatoms (or beads). For example, four water molecules are commonly treated as a single pseudoatom. Such grouping yields two main effects: first, the number of particles is reduced, and second, the time step can be increased, as it is limited by the frequency of the fastest vibrations in the system. However, coarse-grained representation is not restricted to water molecules; Amino Acids, for instance, can also be represented by one or more beads. This dramatically lowers resource requirements, making it feasible to simulate large macromolecular assemblies on timescales up to the micro- and millisecond range. This gain in efficiency, however, inevitably comes at the expense of accuracy compared to full-atom protein descriptions, yielding only semi-quantitative results. Crucial to the success of coarse-grained models is the parameterization of force fields, which must be both accurate and versatile—meaning they are suitable for describing systems of diverse composition and configuration. The larger the beads, the more challenging the parameterization process, because more specific interactions must be effectively captured using a minimal number of parameters and energy terms. This has led to a variety of protein, lipid, and water models representing different compromises between accuracy and universality (see, e.g., Marrink et al. 2004).

Other enhanced sampling Methods based on MD simulations that preserve the atomic representation of the structure include replica exchange molecular dynamics (REMD) and essential dynamics (ED), which are discussed in the following sections. In addition, several non-MD-based methods aimed at predicting protein function are also discussed.

Class="center">9.1.3.1. Replica Exchange Method

The primary goal of most computational modeling studies on biomolecular systems is to calculate macroscopic behavior from microscopic interactions. According to equilibrium Statistical Mechanics, any observable quantity associated with a macroscopic experiment is defined as an ensemble average over all possible states of the system. However, due to the limitations of computing hardware, fully converged sampling of all possible conformational states with their corresponding Boltzmann statistical weights is achievable only for simple systems containing a small number of amino acids (see, e.g., Kubitzki and de Groot 2007). For proteins consisting of hundreds or thousands of amino acids, traditional MD simulations frequently fail to converge, preventing reliable estimation of experimental observables.

Sampling inefficiency is a consequence of the rugged free-energy landscape of the system—a concept introduced by Frauenfelder (Frauenfelder et al. 1991; Frauenfelder and Leeson 1998). Overall, the landscape is generally assumed to be funnel-shaped, with the native states of the system populating the global free-energy minimum (Anfinsen 1973).

Upon closer inspection, the complex multidimensional free-energy landscape is characterized by numerous local minima of nearly equal energy, separated by barriers of varying heights. Each of these minima corresponds to a specific conformational substate, with neighboring minima representing similar Conformations. In terms of this intuitive picture, structural transitions involve overcoming these barriers, and the transition rate depends on the barrier height. In MD simulations at room Temperature, only those barriers that are smaller than or comparable to the thermal energy k_BT can be easily crossed, which corresponds only to minor structural changes, such as side-chain reordering. Consequently, the system spends most of its time trapped in locally stable states (kinetic trapping) rather than exploring diverse conformational states. Such exploration is of immense interest because it relates to biological function, yet it requires the system to overcome high energy barriers. Unfortunately, since MD simulations are mostly limited to the nanosecond timescale, functionally significant conformational transitions are rarely observed.

To tackle this multiple-minima problem, numerous enhanced sampling methods have been proposed (see, e.g., Van Gunsteren and Berendsen 1990; Tai 2004; Adcock and McCammon 2006 and References therein). Among these, generalized-ensemble algorithms have been widely used in recent years (see, e.g., the reviews by Mitsutake et al. 2001; Iba 2001). These algorithms sample an artificial ensemble created by combining or modifying the original ensemble. Algorithms of the second category (e.g., Berg and Neuhaus 1991) primarily alter the initial bell-shaped momentum distribution p(V) of the system by introducing a multicanonical weight factor w(V) such that the resulting distribution becomes flat: p(V)w(V) = const. This flat distribution can then be extensively sampled using MD and Monte Carlo methods, as potential energy barriers are effectively eliminated. Due to these modifications, estimates of physical observables for the canonical ensemble must be obtained via reweighting (Kumar et al. 1992; Chodera et al. 2007). The primary drawback of these algorithms, however, lies in the non-trivial Determination of the various multicanonical weight factors through an iterative Procedure using short test runs. For complex systems, this procedure can be extremely cumbersome, prompting various attempts to improve the convergence of the iterative process (Berg and Celik 1992; Kumar et al. 1996; Smith and Bruce 1996; Hansmann 1997; Bartels and Karplus 1998).

The replica exchange (REX) algorithm, developed as an extension of parallel tempering (Marinari and Parisi 1992), bypasses the weighting factor problem. It belongs to the first category, sampling a generalized ensemble constructed from multiple copies of the original system. Due to its simplicity and ease of Implementation, this algorithm has seen widespread use recently. Most commonly, the standard-temperature formulation of the REX algorithm is employed (Sugita and Okamoto 1999), and the design of the overarching Hamiltonian in this algorithm is receiving increasing attention (Fukunishi et al. 2002; Liu et al. 2005; Sugita et al. 2000; Affentranger et al. 2006; Christen and van Gunsteren 2006; Lyman and Zuckerman 2006).

In replica exchange MD simulations at standard temperature (Sugita and Okamoto 1999), the generalized ensemble is created from M + 1 non-interacting copies, or replicas, of the system spanning a temperature range {T₀,...,T_M} (T_m ≤ T_m+1; m = 0,...,M). This can be implemented, for example, by distributing the calculations across M + 1 nodes on a parallel architecture computer (Fig. 9.6, left). The state of this generalized ensemble is described by a set of states S = {...,s_m,...}, where s_m denotes the state of replica m at temperature T_m. The algorithm consists of two sequential steps: (a) independent simulation of each replica at constant temperature, and (b) swapping replicas S = {...,s_m,...,s_n,...} → S' = {...,s_n,...,s_m,...} according to a Metropolis-like criterion. The acceptance probability for the exchange is given by

where V_m is the potential energy and β_m^-1 = k_BT_m. By alternating steps (a) and (b), the trajectory of the generalized ensemble wanders through temperature space, which in turn leads to diffusion across energy space. This facilitates efficient and statistically rigorous conformational sampling on the system's energy landscape, even in the presence of numerous local minima.

The choice of temperatures strongly influences algorithm performance. Replica temperatures must be selected such that: a) the lowest temperature is sufficiently low for efficient sampling of low-energy states, b) the highest temperature is sufficiently high to overcome energy barriers in the system under study, and c) the acceptance probability P(S→S’) is reasonably high, which requires adequate overlap of the potential energy distributions between neighboring replicas. For large systems with explicitly modeled solvent, this last condition is the most challenging. Simple estimations (Cheng et al. 2005; Fukunishi et al. 2002) show that the largest contribution to the free-energy difference ΔV ~ N_df ΔT comes from the solvent, whose number of degrees of freedom N_{df sol} constitutes the vast majority of the system's total degrees of freedom N_df. Thus, achieving a reasonable acceptance probability can only be accomplished by keeping temperature intervals ΔT = T_m+1 - T_m small (typically a few degrees), which drastically increases the computational requirements for systems with several thousand atoms or more. Despite this strict limitation, REX algorithms have become a standard, widely accepted method for studying peptide folding and denaturation (Zhou et al. 2001; Rao and Caflisch 2003; Garcia and Onuchic 2003; Pitera and Swope 2003; Seibert et al. 2005), structure prediction (Fukunishi et al. 2002; Kokubo and Okamoto 2004), phase transitions (Berg and Neuhaus 1991), and free-energy calculations (Sugita et al. 2000; Lou and Cukier 2006).

Fig. 9.6. Schematic comparison of standard-temperature REX (left) and TEE-REX (right) algorithms for a three-replica simulation. Temperatures are arranged in ascending order, T_i+1 > T_i. Exchange attempts (↔) are undertaken (...) with frequency v_ex. Unlike standard REX, the TEE-REX method excites only the collective subspace {es} (essential subspace) (gray squares), which contains a few collective modes in each replica. The reference replica (T₀, T₀), containing an approximately Boltzmann ensemble, is used for analysis.

To elucidate protein function, another class of enhanced sampling methods that extends beyond traditional MD is also successfully applied. These algorithms build upon the observation that protein fluctuations are generally correlated. The extraction of such collective motion modes and their application to novel sampling algorithms will be the subject of the next two sections.