A cutoff-based method with charge-distribution-data driven pair potentials for efficiently estimating electrostatic interactions in molecular systems.

We introduce a simple cutoff-based method for precise electrostatic energy calculations in the molecular dynamics (MD) simulations of point-particle systems. Our method employs a theoretically derived smooth pair potential function to define electrostatic energy, offering stability and computational efficiency in MD simulations. Instead of imposing specific physical conditions, such as dielectric environments or charge neutrality, we focus on the relationship represented by a single summation formula of charge-weighted pair potentials. This approach allows an accurate energy approximation for each particle, enabling a straightforward error analysis. The resulting particle-dependent pair potential captures the charge distribution information, making it suitable for heterogeneous systems and ensuring an enhanced accuracy through distant information inclusion. Numerical investigations of the Madelung constants of crystalline systems validate the method's accuracy.

Generalized-ensemble method study: A helix-mimetic compound inhibits protein-protein interaction by long-range and short-range intermolecular interactions.

A heterocyclic compound mS-11 is a helix-mimetic designed to inhibit binding of an intrinsic disordered protein neural restrictive silence factor/repressor element 1 silencing factor (NRSF/REST) to a receptor protein mSin3B. We apply a generalized ensemble method, multi-dimensional virtual-system coupled molecular dynamics developed by ourselves recently, to a system consisting of mS-11 and mSin3B, and obtain a thermally equilibrated distribution, which is comprised of the bound and unbound states extensively. The lowest free-energy position of mS-11 coincides with the NRSF/REST position in the experimentally-determined NRSF/REST-mSin3B complex. Importantly, the molecular orientation of mS-11 is ordering in a wide region around mSin3B. The resultant binding scenario is: When mS-11 is distant from the binding site of mSin3B, mS-11 descends the free-energy slope toward the binding site maintaining the molecular orientation to be advantageous for binding. Then, finally a long and flexible hydrophobic sidechain of mS-11 fits into the binding site, which is the lowest-free-energy complex structure inhibiting NRSF/REST binding to mSin3B.

Flexibility and Cell Permeability of Cyclic Ras-Inhibitor Peptides Revealed by the Coupled Nosé-Hoover Equation.

Quantifying the cell permeability of cyclic peptides is crucial for their rational drug design. However, the reasons remain unclear why a minor chemical modification, such as the difference between Ras inhibitors cyclorasin 9A5 and 9A54, can substantially change a peptide's permeability. To address this question, we performed enhanced sampling simulations of these two 11-mer peptides using the coupled Nosé-Hoover equation (cNH) we recently developed. The present cNH simulations realized temperature fluctuations over a wide range (240-600 K) in a dynamic manner, allowing structural samplings that were well validated by nuclear Overhauser effect measurements. The derived structural ensembles were comprehensively analyzed by all-atom structural clustering, mapping the derived clusters onto principal components (PCs) that characterize the cyclic structure, and calculating cluster-dependent geometric and chemical properties. The planar-open conformation was dominant in aqueous solvent, owing to inclusion of the Trp side chain in the main-chain ring, while the compact-closed conformation, which favors cell permeation due to its compactness and high polarity, was also accessible. Conformation-dependent cell permeability was observed in one of the derived PCs, demonstrating that decreased cell permeability in 9A54 is due to the high free energy barrier separating the two conformations. The origin of the change in free energy surface was determined to be loss of flexibility in the modified residues 2-3, resulting from the increased bulkiness of their side chains. The derived molecular mechanism of cell permeability highlights the significance of complete structural dynamics surveys for accelerating drug development with cyclic peptides.

Molecular Interaction Mechanism of a 14-3-3 Protein with a Phosphorylated Peptide Elucidated by Enhanced Conformational Sampling.

Enhanced conformational sampling, a genetic-algorithm-guided multidimensional virtual-system coupled molecular dynamics, can provide equilibrated conformational distributions of a receptor protein and a flexible ligand at room temperature. The distributions provide not only the most stable but also semistable complex structures and propose a ligand-receptor binding process. This method was applied to a system consisting of a receptor protein, 14-3-3ε, and a flexible peptide, phosphorylated myeloid leukemia factor 1 (pMLF1). The results present comprehensive binding pathways of pMLF1 to 14-3-3ε. We identified four thermodynamically stable clusters of MLF1 on the 14-3-3ε surface and free-energy barriers among some clusters. The most stable cluster includes two high-density spots connected by a narrow corridor. When pMLF1 passes the corridor, a salt-bridge relay (switching) related to the phosphorylated residue of pMLF1 occurs. Conformations in one high-density spot are similar to the experimentally determined complex structure. Three-dimensional distributions of residues in the intermolecular interface rationally explain the binding constant changes resulting from the alanine mutation experiment for the residues. We also performed a simulation of nonphosphorylated peptide and 14-3-3ε, which demonstrated that the complex structure was unstable, suggesting that phosphorylation of the peptide is crucially important for binding to 14-3-3ε.

Phosphorylation of an intrinsically disordered region of Ets1 shifts a multi-modal interaction ensemble to an auto-inhibitory state.

Multi-modal interactions are frequently observed in intrinsically disordered regions (IDRs) of proteins upon binding to their partners. In many cases, post-translational modifications in IDRs are accompanied by coupled folding and binding. From both molecular simulations and biochemical experiments with mutational studies, we show that the IDR including a Ser rich region (SRR) of the transcription factor Ets1, just before the DNA-binding core domain, undergoes multi-modal interactions when the SRR is not phosphorylated. In the phosphorylated state, the SRR forms a few specific complex structures with the Ets1 core, covering the recognition helix in the core and drastically reducing the DNA binding affinities as the auto-inhibitory state. The binding kinetics of mutated Ets1 indicates that aromatic residues in the SRR can be substituted with other hydrophobic residues for the interactions with the Ets1 core.

Multidimensional virtual-system coupled canonical molecular dynamics to compute free-energy landscapes of peptide multimer assembly.

An enhanced-sampling method termed multidimensional virtual-system coupled canonical molecular dynamics (mD-VcMD) method is developed. In many cases, generalized-ensemble methods realizing enhanced sampling, for example, adaptive umbrella sampling, apply an effective potential, which is derived from temporarily assumed canonical distribution as a function of one or more arbitrarily defined reaction coordinates. However, it is not straightforward to estimate the appropriate canonical distribution, especially for cases applying multiple reaction coordinates. The current method, mD-VcMD, does not rely on the form of the canonical distribution. Therefore, it is practically useful to explore a high-dimensional reaction-coordinate space. In this article, formulation of mD-VcMD and its evaluation with the simple molecular models consisting of three or four alanine peptides are presented. We confirmed that mD-VcMD efficiently searched 2D and 3D reaction-coordinate spaces defined as interpeptide distances. Direct comparisons with results of long-term canonical MD simulations revealed that mD-VcMD produces correct canonical ensembles. © 2019 Wiley Periodicals, Inc.

Molecular dynamics coupled with a virtual system for effective conformational sampling.

An enhanced conformational sampling method is proposed: virtual-system coupled canonical molecular dynamics (VcMD). Although VcMD enhances sampling along a reaction coordinate, this method is free from estimation of a canonical distribution function along the reaction coordinate. This method introduces a virtual system that does not necessarily obey a physical law. To enhance sampling the virtual system couples with a molecular system to be studied. Resultant snapshots produce a canonical ensemble. This method was applied to a system consisting of two short peptides in an explicit solvent. Conventional molecular dynamics simulation, which is ten times longer than VcMD, was performed along with adaptive umbrella sampling. Free-energy landscapes computed from the three simulations mutually converged well. The VcMD provided quicker association/dissociation motions of peptides than the conventional molecular dynamics did. The VcMD method is applicable to various complicated systems because of its methodological simplicity. © 2018 Wiley Periodicals, Inc.

Multi-dimensional virtual system introduced to enhance canonical sampling.

When an important process of a molecular system occurs via a combination of two or more rare events, which occur almost independently to one another, computational sampling for the important process is difficult. Here, to sample such a process effectively, we developed a new method, named the "multi-dimensional Virtual-system coupled Monte Carlo (multi-dimensional-VcMC)" method, where the system interacts with a virtual system expressed by two or more virtual coordinates. Each virtual coordinate controls sampling along a reaction coordinate. By setting multiple reaction coordinates to be related to the corresponding rare events, sampling of the important process can be enhanced. An advantage of multi-dimensional-VcMC is its simplicity: Namely, the conformation moves widely in the multi-dimensional reaction coordinate space without knowledge of canonical distribution functions of the system. To examine the effectiveness of the algorithm, we introduced a toy model where two molecules (receptor and its ligand) bind and unbind to each other. The receptor has a deep binding pocket, to which the ligand enters for binding. Furthermore, a gate is set at the entrance of the pocket, and the gate is usually closed. Thus, the molecular binding takes place via the two events: ligand approach to the pocket and gate opening. In two-dimensional (2D)-VcMC, the two molecules exhibited repeated binding and unbinding, and an equilibrated distribution was obtained as expected. A conventional canonical simulation, which was 200 times longer than 2D-VcMC, failed in sampling the binding/unbinding effectively. The current method is applicable to various biological systems.

Enhancement of canonical sampling by virtual-state transitions.

A novel method was developed to enhance canonical sampling. A system is divided into virtually introduced sub-states, called "virtual states," which does not exist in reality. The configuration sampling is achieved by a standard canonical sampling method, the Metropolis Monte Carlo method, and confined in a virtual state for a while. In contrast, inter-virtual state motions are controlled by transition probabilities, which can be set arbitrarily. A simple recursive equation was introduced to determine the inter-virtual state transition probabilities, by which the sampling is enhanced considerably. We named this method "virtual-system coupled canonical Monte Carlo (VcMC) sampling." A simple method was proposed to reconstruct a canonical distribution function at a certain temperature from the resultant VcMC sampling data. Two systems, a one-dimensional double-well potential and a three-dimensional ligand-receptor binding/unbinding model, were examined. VcMC produced an accurate canonical distribution much more quickly than a conventional canonical Monte Carlo simulation does.

Enhanced conformational sampling to visualize a free-energy landscape of protein complex formation.

We introduce various, recently developed, generalized ensemble methods, which are useful to sample various molecular configurations emerging in the process of protein-protein or protein-ligand binding. The methods introduced here are those that have been or will be applied to biomolecular binding, where the biomolecules are treated as flexible molecules expressed by an all-atom model in an explicit solvent. Sampling produces an ensemble of conformations (snapshots) that are thermodynamically probable at room temperature. Then, projection of those conformations to an abstract low-dimensional space generates a free-energy landscape. As an example, we show a landscape of homo-dimer formation of an endothelin-1-like molecule computed using a generalized ensemble method. The lowest free-energy cluster at room temperature coincided precisely with the experimentally determined complex structure. Two minor clusters were also found in the landscape, which were largely different from the native complex form. Although those clusters were isolated at room temperature, with rising temperature a pathway emerged linking the lowest and second-lowest free-energy clusters, and a further temperature increment connected all the clusters. This exemplifies that the generalized ensemble method is a powerful tool for computing the free-energy landscape, by which one can discuss the thermodynamic stability of clusters and the temperature dependence of the cluster networks.

Dynamics of the Extended String-Like Interaction of TFIIE with the p62 Subunit of TFIIH.

General transcription factor II E (TFIIE) contains an acid-rich region (residues 378-393) in its α-subunit, comprising 13 acidic and two hydrophobic (Phe387 and Val390) residues. Upon binding to the p62 subunit of TFIIH, the acidic region adopts an extended string-like structure on the basic groove of the pleckstrin homology domain (PHD) of p62, and inserts Phe387 and Val390 into two shallow pockets in the groove. Here, we have examined the dynamics of this interaction by NMR and molecular dynamics (MD) simulations. Although alanine substitution of Phe387 and/or Val390 greatly reduced binding to PHD, the binding mode of the mutants was similar to that of the wild-type, as judged by the chemical-shift changes of the PHD. NMR relaxation dispersion profiles of the interaction exhibited large amplitudes for residues in the C-terminal half-string in the acidic region (Phe387, Glu388, Val390, Ala391, and Asp392), indicating a two-site binding mode: one corresponding to the final complex structure, and one to an off-pathway minor complex. To probe the off-pathway complex structure, an atomically detailed free-energy landscape of the binding mode was computed by all-atom multicanonical MD. The most thermodynamically stable cluster corresponded to the final complex structure. One of the next stable clusters was the off-pathway structure cluster, showing the reversed orientation of the C-terminal half-string on the PHD groove, as compared with the final structure. MD calculations elucidated that the C-terminal half-acidic-string forms encounter complexes mainly around the positive groove region with nearly two different orientations of the string, parallel and antiparallel to the final structure. Interestingly, the most encountered complexes exhibit a parallel-like orientation, suggesting that the string has a tendency to bind around the groove in the proper orientation with the aid of Phe387 and/or Val390 to proceed smoothly to the final complex structure.

Intrinsic disorder accelerates dissociation rather than association.

The intrinsically disordered protein (IDP) has distinct properties both physically and biologically: it often becomes folded when binding to the target and is frequently involved in signal transduction. The physical property seems to be compatible with the biological property where fast association and dissociation between IDP and the target are required. While fast association has been well studied, fueled by the fly-casting mechanism, the dissociation kinetics has received less attention. We here study how the intrinsic disorder affects the dissociation kinetics, as well as the association kinetics, paying attention to the interaction strength at the binding site (i.e., the quality of the "fly lure"). Coarse-grained molecular dynamics simulation of the pKID-KIX system, a well-studied IDP system, shows that the association rate becomes larger as the disorder-inducing flexibility that was imparted to the model is increased, but the acceleration is marginal and turns into deceleration as the quality of the fly lure is worsened. In contrast, the dissociation rate is greatly enhanced as the disorder is increased, indicating that intrinsic disorder serves for rapid signal switching more effectively through dissociation than association. Proteins 2016; 84:1124-1133. © 2016 Wiley Periodicals, Inc.

Variation of free-energy landscape of the p53 C-terminal domain induced by acetylation: Enhanced conformational sampling.

The C-terminal domain (CTD) of tumor suppressor protein p53 is an intrinsically disordered region that binds to various partner proteins, where lysine of CTD is acetylated/nonacetylated and histidine neutralized/non-neutralized. Because of the flexibility of the unbound CTD, a free-energy landscape (FEL) is a useful quantity for determining its statistical properties. We conducted enhanced conformational sampling of CTD in the unbound state via virtual system coupled multicanonical molecular dynamics, in which the lysine was acetylated or nonacetylated and histidine was charged or neutralized. The fragments were expressed by an all-atom model and were immersed in an explicit solvent. The acetylation and charge-neutralization varied FEL greatly, which might be convenient to exert a hub property. The acetylation slightly enhanced alpha-helix structures that are more compact than sheet/loop conformations. The charge-neutralization produced hairpins. Additionally, circular dichroism experiments confirmed the computational results. We propose possible binding mechanisms of CTD to partners by investigating FEL. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.

Virtual-system-coupled adaptive umbrella sampling to compute free-energy landscape for flexible molecular docking.

A novel enhanced conformational sampling method, virtual-system-coupled adaptive umbrella sampling (V-AUS), was proposed to compute 300-K free-energy landscape for flexible molecular docking, where a virtual degrees of freedom was introduced to control the sampling. This degree of freedom interacts with the biomolecular system. V-AUS was applied to complex formation of two disordered amyloid-ß (Aß30-35 ) peptides in a periodic box filled by an explicit solvent. An interpeptide distance was defined as the reaction coordinate, along which sampling was enhanced. A uniform conformational distribution was obtained covering a wide interpeptide distance ranging from the bound to unbound states. The 300-K free-energy landscape was characterized by thermodynamically stable basins of antiparallel and parallel ß-sheet complexes and some other complex forms. Helices were frequently observed, when the two peptides contacted loosely or fluctuated freely without interpeptide contacts. We observed that V-AUS converged to uniform distribution more effectively than conventional AUS sampling did.

Computational Study of Drug Binding Affinity to Influenza A Neuraminidase Using Smooth Reaction Path Generation (SRPG) Method.

Assessment of accurate drug binding affinity to a protein remains a challenge for in silico drug development. In this research, we used the smooth reaction path generation (SRPG) method to calculate binding free energies and determine potential of mean forces (PMFs) along the smoothed dissociation paths of influenza A neuraminidase and its variants with oseltamivir (Tamiflu) and zanamivir (Relenza) inhibitors. With the gained results, we found that the binding free energies of neuraminidase A/H5N1 in WT and two mutants (including H274Y and N294S) with oseltamivir and zanamivir show good agreement with experimental results. Additionally, the thermodynamic origin of the drug resistance of the mutants was also discussed from the PMF profiles.

Multi-scale ensemble modeling of modular proteins with intrinsically disordered linker regions: application to p53.

In eukaryotic proteins, intrinsically disordered regions (IDRs) are ubiquitous and often exist in linker regions that flank the functional domains of modular proteins, regulating their functions. For detailed structural ensemble modeling of IDRs, we propose a multiscale method for IDRs that possess significant long-range order in modular proteins and apply it to the eukaryotic transcription factor p53 as an example. First, we performed all-atom (AA) molecular dynamics (MD) simulations of the explicitly solvated p53 linker region, without experimental restraint terms, finding fractional long-range contacts within the linker. Second, we fed this AA MD ensemble into a coarse-grained (CG) model, finding an optimal set of contact potentials. The optimized CG MD simulations reproduced the contact probability map from the AA MD simulations. Finally, we performed the CG MD simulation of the tetrameric p53 fragments including the core domains, the linker, and the tetramerization domain. Using the obtained ensemble, we theoretically calculated the small angle x-ray scattering (SAXS) profile of this fragment. The obtained SAXS profile agrees well with the experiment. We also found that the long-range contacts in the p53 linker region are required to reproduce the experimental SAXS profile. The developed framework in which we calculate the long-range contact probability map from the AA MD simulation and incorporate it to the CG model can be applied to broad range of IDRs.

High-resolution modeling of antibody structures by a combination of bioinformatics, expert knowledge, and molecular simulations.

In the second antibody modeling assessment, we used a semiautomated template-based structure modeling approach for 11 blinded antibody variable region (Fv) targets. The structural modeling method involved several steps, including template selection for framework and canonical structures of complementary determining regions (CDRs), homology modeling, energy minimization, and expert inspection. The submitted models for Fv modeling in Stage 1 had the lowest average backbone root mean square deviation (RMSD) (1.06 Å). Comparison to crystal structures showed the most accurate Fv models were generated for 4 out of 11 targets. We found that the successful modeling in Stage 1 mainly was due to expert-guided template selection for CDRs, especially for CDR-H3, based on our previously proposed empirical method (H3-rules) and the use of position specific scoring matrix-based scoring. Loop refinement using fragment assembly and multicanonical molecular dynamics (McMD) was applied to CDR-H3 loop modeling in Stage 2. Fragment assembly and McMD produced putative structural ensembles with low free energy values that were scored based on the OSCAR all-atom force field and conformation density in principal component analysis space, respectively, as well as the degree of consensus between the two sampling methods. The quality of 8 out of 10 targets improved as compared with Stage 1. For 4 out of 10 Stage-2 targets, our method generated top-scoring models with RMSD values of less than 1 Å. In this article, we discuss the strengths and weaknesses of our approach as well as possible directions for improvement to generate better predictions in the future.

Free-energy landscape of intrinsically disordered proteins investigated by all-atom multicanonical molecular dynamics.

We introduce computational studies on intrinsically disordered proteins (IDPs). Especially, we present our multicanonical molecular dynamics (McMD) simulations of two IDP-partner systems: NRSF-mSin3 and pKID-KIX. McMD is one of enhanced conformational sampling methods useful for conformational sampling of biomolecular systems. IDP adopts a specific tertiary structure upon binding to its partner molecule, although it is unstructured in the unbound state (i.e. the free state). This IDP-specific property is called "coupled folding and binding". The McMD simulation treats the biomolecules with an all-atom model immersed in an explicit solvent. In the initial configuration of simulation, IDP and its partner molecules are set to be distant from each other, and the IDP conformation is disordered. The computationally obtained free-energy landscape for coupled folding and binding has shown that native- and non-native-complex clusters distribute complicatedly in the conformational space. The all-atom simulation suggests that both of induced-folding and population-selection are coupled complicatedly in the coupled folding and binding. Further analyses have exemplified that the conformational fluctuations (dynamical flexibility) in the bound and unbound states are essentially important to characterize IDP functioning.

Binding Mechanism of Riboswitch to Natural Ligand Elucidated by McMD-Based Dynamic Docking Simulations.

Flavin mononucleotide riboswitches are common among many pathogenic bacteria and are therefore considered to be an attractive target for antibiotics development. The riboswitch binds riboflavin (RBF, also known as vitamin B2), and although an experimental structure of their complex has been solved with the ligand bound deep inside the RNA molecule in a seemingly unreachable state, the binding mechanism between these molecules is not yet known. We have therefore used our Multicanonical Molecular Dynamics (McMD)-based dynamic docking protocol to analyze their binding mechanism by simulating the binding process between the riboswitch aptamer domain and the RBF, starting from the apo state of the riboswitch. Here, the refinement stage was crucial to identify the native binding configuration, as several other binding configurations were also found by McMD-based docking simulations. RBF initially binds the interface between P4 and P6 including U61 and G62, which forms a gateway where the ligand lingers until this gateway opens sufficiently to allow the ligand to pass through and slip into the hidden binding site including A48, A49, and A85.

A virtual-system coupled multicanonical molecular dynamics simulation: principles and applications to free-energy landscape of protein-protein interaction with an all-atom model in explicit solvent.

We propose a novel generalized ensemble method, a virtual-system coupled multicanonical molecular dynamics (V-McMD), to enhance conformational sampling of biomolecules expressed by an all-atom model in an explicit solvent. In this method, a virtual system, of which physical quantities can be set arbitrarily, is coupled with the biomolecular system, which is the target to be studied. This method was applied to a system of an Endothelin-1 derivative, KR-CSH-ET1, known to form an antisymmetric homodimer at room temperature. V-McMD was performed starting from a configuration in which two KR-CSH-ET1 molecules were mutually distant in an explicit solvent. The lowest free-energy state (the most thermally stable state) at room temperature coincides with the experimentally determined native complex structure. This state was separated to other non-native minor clusters by a free-energy barrier, although the barrier disappeared with elevated temperature. V-McMD produced a canonical ensemble faster than a conventional McMD method.