GENESIS 2.1: High-Performance Molecular Dynamics Software for Enhanced Sampling and Free-Energy Calculations for Atomistic, Coarse-Grained, and Quantum Mechanics/Molecular Mechanics Models.

GENeralized-Ensemble SImulation System (GENESIS) is a molecular dynamics (MD) software developed to simulate the conformational dynamics of a single biomolecule, as well as molecular interactions in large biomolecular assemblies and between multiple biomolecules in cellular environments. To achieve the latter purpose, the earlier versions of GENESIS emphasized high performance in atomistic MD simulations on massively parallel supercomputers, with or without graphics processing units (GPUs). Here, we implemented multiscale MD simulations that include atomistic, coarse-grained, and hybrid quantum mechanics/molecular mechanics (QM/MM) calculations. They demonstrate high performance and are integrated with enhanced conformational sampling algorithms and free-energy calculations without using external programs except for the QM programs. In this article, we review new functions, molecular models, and other essential features in GENESIS version 2.1 and discuss ongoing developments for future releases.

Modified Protein-Water Interactions in CHARMM36m for Thermodynamics and Kinetics of Proteins in Dilute and Crowded Solutions.

Proper balance between protein-protein and protein-water interactions is vital for atomistic molecular dynamics (MD) simulations of globular proteins as well as intrinsically disordered proteins (IDPs). The overestimation of protein-protein interactions tends to make IDPs more compact than those in experiments. Likewise, multiple proteins in crowded solutions are aggregated with each other too strongly. To optimize the balance, Lennard-Jones (LJ) interactions between protein and water are often increased about 10% (with a scaling parameter, λ = 1.1) from the existing force fields. Here, we explore the optimal scaling parameter of protein-water LJ interactions for CHARMM36m in conjunction with the modified TIP3P water model, by performing enhanced sampling MD simulations of several peptides in dilute solutions and conventional MD simulations of globular proteins in dilute and crowded solutions. In our simulations, 10% increase of protein-water LJ interaction for the CHARMM36m cannot maintain stability of a small helical peptide, (AAQAA)3 in a dilute solution and only a small modification of protein-water LJ interaction up to the 3% increase (λ = 1.03) is allowed. The modified protein-water interactions are applicable to other peptides and globular proteins in dilute solutions without changing thermodynamic properties from the original CHARMM36m. However, it has a great impact on the diffusive properties of proteins in crowded solutions, avoiding the formation of too sticky protein-protein interactions.

Modified Hamiltonian in FEP Calculations for Reducing the Computational Cost of Electrostatic Interactions.

The free-energy perturbation (FEP) method predicts relative and absolute free-energy changes of biomolecules in solvation and binding with other molecules. FEP is, therefore, one of the most essential tools in in silico drug design. In conventional FEP, to smoothly connect two thermodynamic states, the potential energy is modified as a linear combination of the end-state potential energies by introducing scaling factors. When the particle mesh Ewald is used for electrostatic calculations, conventional FEP requires two reciprocal-space calculations per time step, which largely decreases the computational performance. To overcome this problem, we propose a new FEP scheme by introducing a modified Hamiltonian instead of interpolation of the end-state potential energies. The scheme introduces nonuniform scaling into the electrostatic potential as used in Replica Exchange with Solute Tempering 2 (REST2) and does not require additional reciprocal-space calculations. We tested this modified Hamiltonian in FEP calculations in several biomolecular systems. In all cases, the calculated free-energy changes with the current scheme are in good agreement with those from conventional FEP. The modified Hamiltonian in FEP greatly improves the computational performance, which is particularly marked for large biomolecular systems whose reciprocal-space calculations are the major bottleneck of total computational time.

Use of multistate Bennett acceptance ratio method for free-energy calculations from enhanced sampling and free-energy perturbation.

Multistate Bennett acceptance ratio (MBAR) works as a method to analyze molecular dynamics (MD) simulation data after the simulations have been finished. It is widely used to estimate free-energy changes between different states and averaged properties at the states of interest. MBAR allows us to treat a wide range of states from those at different temperature/pressure to those with different model parameters. Due to the broad applicability, the MBAR equations are rather difficult to apply for free-energy calculations using different types of MD simulations including enhanced conformational sampling methods and free-energy perturbation. In this review, we first summarize the basic theory of the MBAR equations and categorize the representative usages into the following four: (i) perturbation, (ii) scaling, (iii) accumulation, and (iv) full potential energy. For each, we explain how to prepare input data using MD simulation trajectories for solving the MBAR equations. MBAR is also useful to estimate reliable free-energy differences using MD trajectories based on a semi-empirical quantum mechanics/molecular mechanics (QM/MM) model and ab initio QM/MM energy calculations on the MD snapshots. We also explain how to use the MBAR software in the GENESIS package, which we call mbar_analysis, for the four representative cases. The proposed estimations of free-energy changes and thermodynamic averages are effective and useful for various biomolecular systems.

Optimized Hydrogen Mass Repartitioning Scheme Combined with Accurate Temperature/Pressure Evaluations for Thermodynamic and Kinetic Properties of Biological Systems.

In recent years, molecular dynamics (MD) simulations with larger time steps have been performed with the hydrogen-mass-repartitioning (HMR) scheme, where the mass of each hydrogen atom is increased to reduce high-frequency motion while the mass of a non-hydrogen atom bonded to a hydrogen atom is decreased to keep the total molecular mass unchanged. Here, we optimize the scaling factors in HMR and combine them with previously developed accurate temperature/pressure evaluations. The heterogeneous HMR scaling factors are useful to avoid the structural instability of amino acid residues having a five- or six-membered ring in MD simulations with larger time steps. It also reproduces kinetic properties, namely translational and rotational diffusions, if the HMR scaling factors are applied to only solute molecules. The integration scheme is tested for biological systems that include soluble/membrane proteins and lipid bilayers for about 200 µs MD simulations in total and give consistent results in MD simulations with both a small time step of 2.0 fs and a large, multiple time step integration with time steps of 3.5 fs (for fast motions) and 7.0 fs (for slower motions). We also confirm that the multiple time step integration scheme used in this study provides more accurate energy conservations than the RESPA/C1 and is comparable to the RESPA/C2 in NAMD. In summary, the current integration scheme combining the optimized HMR with accurate temperature/pressure evaluations can provide stable and reliable MD trajectories with a larger time step, which are computationally more than 2-fold efficient compared to the conventional methods.

Reduced efficacy of a Src kinase inhibitor in crowded protein solution.

The inside of a cell is highly crowded with proteins and other biomolecules. How proteins express their specific functions together with many off-target proteins in crowded cellular environments is largely unknown. Here, we investigate an inhibitor binding with c-Src kinase using atomistic molecular dynamics (MD) simulations in dilute as well as crowded protein solution. The populations of the inhibitor, 4-amino-5-(4-methylphenyl)-7-(t-butyl)pyrazolo[3,4-d]pyrimidine (PP1), in bulk solution and on the surface of c-Src kinase are reduced as the concentration of crowder bovine serum albumins (BSAs) increases. This observation is consistent with the reduced PP1 inhibitor efficacy in experimental c-Src kinase assays in addition with BSAs. The crowded environment changes the major binding pathway of PP1 toward c-Src kinase compared to that in dilute solution. This change is explained based on the population shift mechanism of local conformations near the inhibitor binding site in c-Src kinase.

Unraveling the Coupling between Conformational Changes and Ligand Binding in Ribose Binding Protein Using Multiscale Molecular Dynamics and Free-Energy Calculations.

Conformational changes of proteins upon ligand binding are usually explained in terms of several mechanisms including the induced fit, conformational selection, or their mixtures. Due to the slow time scales, conventional molecular dynamics (cMD) simulations based on the atomistic models cannot easily simulate the open-to-closed conformational transition in proteins. In our previous study, we have developed an enhanced sampling scheme (generalized replica exchange with solute tempering selected surface charged residues: gREST_SSCR) for multidomain proteins and applied it to ligand-mediated conformational changes in the G134R mutant of ribose-binding protein (RBPG134R) in solution. The free-energy landscape (FEL) of RBPG134R in the presence of a ribose at the binding site included the open and closed states and two intermediates, open-like and closed-like forms. Only the open and open-like forms existed in the FEL without a ribose. In the current study, the coupling between the conformational changes and ligand binding is further investigated using coarse-grained MD, multiple atomistic cMD, and free-energy calculations. The ribose is easily dissociated from the binding site of wild-type RBP and RBPG134R in the cMD simulations starting from the open and open-like forms. In contrast, it is stable at the binding site in the simulations from the closed and closed-like forms. The free-energy calculations provide the binding affinities of different structures, supporting the results of cMD simulations. Importantly, cMD simulations from the closed-like structures reveal transitions toward the closed one in the presence of a bound ribose. On the basis of the computational results, we propose a molecular mechanism in which conformational selection and induced fit happen in the first and second halves of the open-to-closed transition in RBP, respectively.

CHARMM-GUI Free Energy Calculator for Absolute and Relative Ligand Solvation and Binding Free Energy Simulations.

Alchemical free energy simulations have long been utilized to predict free energy changes for binding affinity and solubility of small molecules. However, while the theoretical foundation of these methods is well established, seamlessly handling many of the practical aspects regarding the preparation of the different thermodynamic end states of complex molecular systems and the numerous processing scripts often remains a burden for successful applications. In this work, we present CHARMM-GUI Free Energy Calculator (http://www.charmm-gui.org/input/fec) that provides various alchemical free energy perturbation molecular dynamics (FEP/MD) systems with input and post-processing scripts for NAMD and GENESIS. Four submodules are available: Absolute Ligand Binder (for absolute ligand binding FEP/MD), Relative Ligand Binder (for relative ligand binding FEP/MD), Absolute Ligand Solvator (for absolute ligand solvation FEP/MD), and Relative Ligand Solvator (for relative ligand solvation FEP/MD). Each module is designed to build multiple systems of a set of selected ligands at once for high-throughput FEP/MD simulations. The capability of Free Energy Calculator is illustrated by absolute and relative solvation FEP/MD of a set of ligands and absolute and relative binding FEP/MD of a set of ligands for T4-lysozyme in solution and the adenosine A2A receptor in a membrane. The calculated free energy values are overall consistent with the experimental and published free energy results (within ∼1 kcal/mol). We hope that Free Energy Calculator is useful to carry out high-throughput FEP/MD simulations in the field of biomolecular sciences and drug discovery.

Prediction of Protein-Ligand Binding Pose and Affinity Using the gREST+FEP Method.

The accurate prediction of protein-ligand binding affinity is a central challenge in computational chemistry and in-silico drug discovery. The free energy perturbation (FEP) method based on molecular dynamics (MD) simulation provides reasonably accurate results only if a reliable structure is available via high-resolution X-ray crystallography. To overcome the limitation, we propose a sequential prediction protocol using generalized replica exchange with solute tempering (gREST) and FEP. At first, ligand binding poses are predicted using gREST, which weakens protein-ligand interactions at high temperatures to sample multiple binding poses. To avoid ligand dissociation at high temperatures, a flat-bottom restraint potential centered on the binding site is applied in the simulation. The binding affinity of the most reliable pose is then calculated using FEP. The protocol is applied to the bindings of ten ligands to FK506 binding proteins (FKBP), showing the excellent agreement between the calculated and experimental binding affinities. The present protocol, which is referred to as the gREST+FEP method, would help to predict the binding affinities without high-resolution structural information on the ligand-bound state.

Mosaic Cooperativity in Slow Polypeptide Topological Isomerization Revealed by Residue-Specific NMR Thermodynamic Analysis.

Slow polypeptide conformational changes on time scales of >1 s are generally assumed to be highly cooperative two-state transitions, reflecting the high energy barrier. However, few experimental characterizations have tested the validity of this assumption. We performed residue-specific NMR thermodynamic analysis of the 27-residue lantibiotic peptide, nukacin ISK-1, to characterize the isomerization between two topological states on the second time scale. Unexpectedly, the thermal transition behaviors were distinct among peptide regions, indicating that the topological isomerization process is a mosaic of different degrees of cooperativity. The conformational change path between the two NMR structures was deduced by a targeted molecular dynamics simulation. The unique side-chain threading motions through the monosulfide rings are the structural basis of the high energy barrier, and the nonlocal interactions in the hydrophobic core are the structural basis of the cooperativity. Taken together, we provide an energetic description of the topological isomerization of nukacin ISK-1.

Structural Mechanisms Underlying Activity Changes in an AMPA-type Glutamate Receptor Induced by Substitutions in Its Ligand-Binding Domain.

α-Amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA)-type glutamate receptors produce postsynaptic current by transmitting an agonist-induced structural change in the ligand-binding domain (LBD) to the transmembrane channel. Receptors carrying T686S/A substitutions in their LBDs produce weaker glutamate-evoked currents than wild-type (WT) receptors. However, the substitutions induce little differences in the crystal structures of their LBDs. To understand the structural mechanism underlying reduced activities of these AMPAR variants, we analyzed the structural dynamics of WT, T686S, and T686A variants of LBD using nuclear magnetic resonance. The HD exchange studies of the LBDs showed that the kinetic step where the ligand-binding cleft closes was changed by the substitutions, and the substitution-induced population shift from cleft-closed to cleft-open structures is responsible for the reduced activities of the variants. The chemical shift analyses revealed another structural equilibrium between cleft-locked and cleft-partially-open conformations. The substitution-induced population shift in this equilibrium may be related to slower desensitization observed for these variants.

Replica-Exchange Umbrella Sampling Combined with Gaussian Accelerated Molecular Dynamics for Free-Energy Calculation of Biomolecules.

We have developed an enhanced conformational sampling method combining replica-exchange umbrella sampling (REUS) with Gaussian accelerated molecular dynamics (GaMD). REUS enhances the sampling along predefined reaction coordinates, while GaMD accelerates the conformational dynamics by adding a boost potential to the system energy. The method, which we call GaREUS (Gaussian accelerated replica-exchange umbrella sampling), enhances the sampling more efficiently than REUS or GaMD, while the computational resource for GaREUS is the same as that required for REUS. The two-step reweighting procedure using the multistate Bennett acceptance ratio method and the cumulant expansion for the exponential average is applied to the simulation trajectories for obtaining the unbiased free-energy landscapes. We apply GaREUS to the calculations of free-energy landscapes for three different cases: conformational equilibria of N-glycan, folding of chignolin, and conformational change of adenyl kinase. We show that GaREUS speeds up the convergences of free-energy calculations using the same amount of computational resources as REUS. The free-energy landscapes reweighted from the trajectories of GaREUS agree with previously reported ones. GaREUS is applicable to free-energy calculations of various biomolecular dynamics and functions with reasonable computational costs.

Encounter complexes and hidden poses of kinase-inhibitor binding on the free-energy landscape.

Modern drug discovery increasingly focuses on the drug-target binding kinetics which depend on drug (un)binding pathways. The conventional molecular dynamics simulation can observe only a few binding events even using the fastest supercomputer. Here, we develop 2D gREST/REUS simulation with enhanced flexibility of the ligand and the protein binding site. Simulation (43 µs in total) applied to an inhibitor binding to c-Src kinase covers 100 binding and unbinding events. On the statistically converged free-energy landscapes, we succeed in predicting the X-ray binding structure, including water positions. Furthermore, we characterize hidden semibound poses and transient encounter complexes on the free-energy landscapes. Regulatory residues distant from the catalytic core are responsible for the initial inhibitor uptake and regulation of subsequent bindings, which was unresolved by experiments. Stabilizing/blocking of either the semibound poses or the encounter complexes can be an effective strategy to optimize drug-target residence time.

Replica-Exchange Methods for Biomolecular Simulations.

In this study, a replica-exchange method was developed to overcome conformational sampling difficulties in computer simulations of spin glass or other systems with rugged free-energy landscapes. This method was then applied to the protein-folding problem in combination with molecular dynamics (MD) simulation. Owing to its simplicity and sampling efficiency, the replica-exchange method has been applied to many other biological problems and has been continuously improved. The method has often been combined with other sampling techniques, such as umbrella sampling, free-energy perturbation, metadynamics, and Gaussian accelerated MD (GaMD). In this chapter, we first summarize the original replica-exchange molecular dynamics (REMD) method and discuss how new algorithms related to the original method are implemented to add new features. Heterogeneous and flexible structures of an N-glycan in a solution are simulated as an example of applications by REMD, replica exchange with solute tempering, and GaMD. The sampling efficiency of these methods on the N-glycan system and the convergence of the free-energy changes are compared. REMD simulation protocols and trajectory analysis using the GENESIS software are provided to facilitate the practical use of advanced simulation methods.

De Novo Prediction of Binders and Nonbinders for T4 Lysozyme by gREST Simulations.

Molecular recognition underpins all specific protein-ligand interactions and is essential for biomolecular functions. The prediction of canonical binding poses and distinguishing binders from nonbinders are much sought after goals. Here, we apply the generalized replica exchange with solute tempering method, gREST, combined with a flat-bottom potential to evaluate binder and nonbinder interactions with a T4 lysozyme Leu99Ala mutant. The buried hydrophobic cavity and possibility of coupled conformational changes in this protein make binding predictions difficult. The present gREST simulations, enabling enhanced flexibilities of the ligand and protein residues near the binding site, sample bindings in multiple poses, and correct portrayal of X-ray structures. The free-energy profiles of binders (benzene, ethylbenzene, and n-hexylbenzene) are distinct from those of nonbinders (phenol and benzaldehyde). Bindings of the two larger molecules seem to be associated with a structural change toward an excited conformation of the protein, which agrees with experimental findings. The protocol is generally applicable to various proteins having buried cavities with limited access for ligands with different shapes, sizes, and chemical properties.

Population Shift Mechanism for Partial Agonism of AMPA Receptor.

α-amino-3-hydroxy-5-methyl-4-isoaxazolepropionic acid (AMPA) ionotropic glutamate receptors mediate fast excitatory neurotransmission in the central nervous system, and their dysfunction is associated with neurological diseases. Glutamate binding to ligand-binding domains (LBDs) of AMPA receptors induces channel opening in the transmembrane domains of the receptors. The T686A mutation reduces glutamate efficacy so that the glutamate behaves as a partial agonist. The crystal structures of wild-type and mutant LBDs are very similar and cannot account for the observed behavior. To elucidate the molecular mechanism inducing partial agonism of the T686A mutant, we computed the free-energy landscapes governing GluA2 LBD closure using replica-exchange umbrella sampling simulations. A semiclosed state, not observed in crystal structures, appears in the mutant during simulation. In this state, the LBD cleft opens slightly because of breaking of interlobe hydrogen bonds, reducing the efficiency of channel opening. The energy difference between the LBD closed and semiclosed states is small, and transitions between the two states would occur by thermal fluctuations. Evidently, glutamate binding to the T686A mutant induces a population shift from a closed to a semiclosed state, explaining the partial agonism in the AMPA receptor.

Water based on a molecular model behaves like a hard-sphere solvent for a nonpolar solute when the reference interaction site model and related theories are employed.

For neutral hard-sphere solutes, we compare the reduced density profile of water around a solute g(r), solvation free energy µ, energy U, and entropy S under the isochoric condition predicted by the two theories: dielectrically consistent reference interaction site model (DRISM) and angle-dependent integral equation (ADIE) theories. A molecular model for water pertinent to each theory is adopted. The hypernetted-chain (HNC) closure is employed in the ADIE theory, and the HNC and Kovalenko-Hirata (K-H) closures are tested in the DRISM theory. We also calculate g(r), U, S, and µ of the same solute in a hard-sphere solvent whose molecular diameter and number density are set at those of water, in which case the radial-symmetric integral equation (RSIE) theory is employed. The dependences of µ, U, and S on the excluded volume and solvent-accessible surface area are analyzed using the morphometric approach (MA). The results from the ADIE theory are in by far better agreement with those from computer simulations available for g(r), U, and µ. For the DRISM theory, g(r) in the vicinity of the solute is quite high and becomes progressively higher as the solute diameter d U increases. By contrast, for the ADIE theory, it is much lower and becomes further lower as d U increases. Due to unphysically positive U and significantly larger |S|, µ from the DRISM theory becomes too high. It is interesting that µ, U, and S from the K-H closure are worse than those from the HNC closure. Overall, the results from the DRISM theory with a molecular model for water are quite similar to those from the RSIE theory with the hard-sphere solvent. Based on the results of the MA analysis, we comparatively discuss the different theoretical methods for cases where they are applied to studies on the solvation of a protein.

Statistical Thermodynamics for Actin-Myosin Binding: The Crucial Importance of Hydration Effects.

Actomyosin is an important molecular motor, and the binding of actin and myosin is an essential research target in biophysics. Nevertheless, the physical factors driving or opposing the binding are still unclear. Here, we investigate the role of water in actin-myosin binding using the most reliable statistical-mechanical method currently available for assessing biomolecules immersed in water. This method is characterized as follows: water is treated not as a dielectric continuum but as an ensemble of molecules; the polyatomic structures of proteins are taken into consideration; and the binding free energy is decomposed into physically insightful entropic and energetic components by accounting for the hydration effect to its full extent. We find that the actin-myosin binding brings large gains of electrostatic and Lennard-Jones attractive interactions. However, these gains are accompanied by even larger losses of actin-water and myosin-water electrostatic and LJ attractive interactions. Although roughly half of the energy increase due to the losses is cancelled out by the energy decrease arising from structural reorganization of the water released upon binding, the remaining energy increase is still larger than the energy decrease brought by the gains mentioned above. Hence, the net change in system energy is positive, which opposes binding. Importantly, the binding is driven by a large gain of configurational entropy of water, which surpasses the positive change in system energy and the conformational entropy loss occurring for actin and myosin. The principal physical origin of the large water-entropy gain is as follows: the actin-myosin interface is closely packed with the achievement of high shape complementarity on the atomic level, leading to a large increase in the total volume available to the translational displacement of water molecules in the system and a resultant reduction of water crowding (i.e., entropic correlations among water molecules).

A highly efficient hybrid method for calculating the hydration free energy of a protein.

We develop a new method for calculating the hydration free energy (HFE) of a protein with any net charge. The polar part of the energetic component in the HFE is expressed as a linear combination of four geometric measures (GMs) of the protein structure and the generalized Born (GB) energy plus a constant. The other constituents in the HFE are expressed as linear combinations of the four GMs. The coefficients (including the constant) in the linear combinations are determined using the three-dimensional reference interaction site model (3D-RISM) theory applied to sufficiently many protein structures. Once the coefficients are determined, the HFE and its constituents of any other protein structure are obtained simply by calculating the four GMs and GB energy. Our method and the 3D-RISM theory give perfectly correlated results. Nevertheless, the computation time required in our method is over four orders of magnitude shorter.