Molecular dynamics scheme for precise estimation of electrostatic interaction via zero-dipole summation principle

Analytic Gradients for the Electrostatic Embedding QM/MM Model in Periodic Boundary Conditions Using Particle-Mesh Ewald Sums and Electrostatic Potential Fitted Charge Operators

Long-range electrostatic effects are fundamental for describing chemical reactivity in the condensed phase. Here, we present the methodology of an efficient quantum mechanical/molecular mechanical (QM/MM) model in periodic boundary conditions (PBC) compatible with QM/MM boundaries at chemical bonds. The method combines electrostatic potential fitted charge operators and electrostatic potentials derived from the smooth particle-mesh Ewald (PME) sum approach. The total energy and its analytic first derivatives with respect to QM, MM, and lattice vectors allow QM/MM molecular dynamics (MD) in the most common thermodynamic ensembles. We demonstrate the robustness of the method by performing a QM/MM MD equilibration of methanol in water. We simulate the cis/trans isomerization free-energy profiles in water of proline amino acid and a proline-containing oligopeptide, showing a correct description of the reaction barrier. Our PBC-compatible QM/MM model can efficiently be used to study the chemical reactivity in the condensed phase and enzymatic catalysis.

Binding Mechanism between Platelet Glycoprotein and Cyclic Peptide Elucidated by McMD-Based Dynamic Docking

The cyclic peptide OS1 (amino acid sequence: CTERMALHNLC), which has a disulfide bond between both termini cysteine residues, inhibits complex formation between the platelet glycoprotein Ibα (GPIbα) and the von Willebrand factor (vWF) by forming a complex with GPIbα. To study the binding mechanism between GPIbα and OS1 and, therefore, the inhibition mechanism of the protein-protein GPIbα-vWF complex, we have applied our multicanonical molecular dynamics (McMD)-based dynamic docking protocol starting from the unbound state of the peptide. Our simulations have reproduced the experimental complex structure, although the top-ranking structure was an intermediary one, where the peptide was bound in the same location as in the experimental structure; however, the β-switch of GPIbα attained a different conformation. Our analysis showed that subsequent refolding of the β-switch results in a more stable binding configuration, although the transition to the native configuration appears to take some time, during which OS1 could dissociate. Our results show that conformational changes in the β-switch are crucial for successful binding of OS1. Furthermore, we identified several allosteric binding sites of GPIbα that might also interfere with vWF binding, and optimization of the peptide to target these allosteric sites might lead to a more effective inhibitor, as these are not dependent on the β-switch conformation.

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.

Molecular dynamics simulation analysis of structural dynamic cross correlation induced by odorant hydrogen-bonding in mouse eugenol ol- factory receptor

Structural fluctuations and dynamic cross-correlations in the mouse eugenol olfactory receptor (Olfr73) were studied by molecular dynamics (MD) simulation to characterize the dynamic response of the protein upon ligand binding. The initial structure was generated by the artificial intelligence tool AlphaFold2 due to the current lack of experimental data. We focused on the hydrogen (H) bond of the odorant eugenol to Ser113, Asn207, and Tyr260 of the receptor protein, the importance of which has been suggested by previous experimental studies. The H-bond was not observed in docking simulations, but in subsequent MD simulations the H-bond to Ser113 was formed in 2-4 ns. The lifetime of the H-bond was in the range of 1-20 ns. On the trajectory with the most stable (20 ns) H-bond, the structural fluctuation of the α-carbon atoms of the receptor main chain was studied by calculating the root mean square fluctuations, the dynamic cross-correlation map, and the time-dependent dynamic cross-correlation. The analysis suggested a correlation transfer pathway Ser113 → Phe182 → (Leu259 or Tyr260) → Tyr291 induced by the ligand binding with a time scale of 4-6 ns.

Enhanced Coarse-Grained Molecular Dynamics Simulation with a Smoothed Hybrid Potential Using a Neural Network Model

In all-atom (AA) molecular dynamics (MD) simulations, the rugged energy profile of the force field makes it challenging to reproduce spontaneous structural changes in biomolecules within a reasonable calculation time. Existing coarse-grained (CG) models, in which the energy profile is set to a global minimum around the initial structure, are unsuitable to explore the structural dynamics between metastable states far away from the initial structure without any bias. In this study, we developed a new hybrid potential composed of an artificial intelligence (AI) potential and minimal CG potential related to the statistical bond length and excluded volume interactions to accelerate the transition dynamics while maintaining the protein character. The AI potential is trained by energy matching using a diverse structural ensemble sampled via multicanonical (Mc) MD simulation and the corresponding AA force field energy, profile of which is smoothed by energy minimization. By applying the new methodology to chignolin and TrpCage, we showed that the AI potential can predict the AA energy with significantly high accuracy, as indicated by a correlation coefficient (R-value) between the true and predicted energies exceeding 0.89. In addition, we successfully demonstrated that CGMD simulation based on the smoothed hybrid potential can significantly enhance the transition dynamics between various metastable states while preserving protein properties compared to those obtained with conventional CGMD and AAMD.

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.

Elucidation of binding mechanism, affinity, and complex structure between mWT1 tumor-associated antigen peptide and HLA-A*24:02

We have applied our advanced computational and experimental methodologies to investigate the complex structure and binding mechanism of a modified Wilms' Tumor 1 (mWT1) protein epitope to the understudied Asian-dominant allele HLA-A*24:02 (HLA-A24) in aqueous solution. We have applied our developed multicanonical molecular dynamics (McMD)-based dynamic docking method to analyze the binding pathway and mechanism, which we verified by comparing the highest probability structures from simulation with our experimentally solved x-ray crystal structure. Subsequent path sampling MD simulations elucidated the atomic details of the binding process and indicated that first an encounter complex is formed between the N-terminal's positive charge of the 9-residue mWT1 fragment peptide and a cluster of negative residues on the surface of HLA-A24, with the major histocompatibility complex (MHC) molecule preferring a predominantly closed conformation. The peptide first binds to this closed MHC conformation, forming an encounter complex, after which the binding site opens due to increased entropy of the binding site, allowing the peptide to bind to form the native complex structure. Further sequence and structure analyses also suggest that although the peptide loading complex would help with stabilizing the MHC molecule, the binding depends in a large part on the intrinsic affinity between the MHC molecule and the antigen peptide. Finally, our computational tools and analyses can be of great benefit to study the binding mechanism of different MHC types to their antigens, where it could also be useful in the development of higher affinity variant peptides and for personalized medicine.

Mutual induced-fit mechanism drives binding between intrinsically disordered Bim and cryptic binding site of Bcl-xL

The intrinsically disordered region (IDR) of Bim binds to the flexible cryptic site of Bcl-xL, a pro-survival protein involved in cancer progression that plays an important role in initiating apoptosis. However, their binding mechanism has not yet been elucidated. We have applied our dynamic docking protocol, which correctly reproduced both the IDR properties of Bim and the native bound configuration, as well as suggesting other stable/meta-stable binding configurations and revealed the binding pathway. Although the cryptic site of Bcl-xL is predominantly in a closed conformation, initial binding of Bim in an encounter configuration leads to mutual induced-fit binding, where both molecules adapt to each other; Bcl-xL transitions to an open state as Bim folds from a disordered to an α-helical conformation while the two molecules bind each other. Finally, our data provides new avenues to develop novel drugs by targeting newly discovered stable conformations of Bcl-xL.

Thermal Stability Estimation of Single Domain Antibodies Using Molecular Dynamics Simulations

In this chapter, we describe a protocol to estimate the thermal stability of single domain antibodies (sdAbs) using molecular dynamics (MD) simulations. This method measures the Q-value, the fraction of the native contacts, along the trajectory of high-temperature MD simulations starting from the experimental X-ray structure. We show a good correlation between the Q-value and the experimental melting temperature (T_m) in seven sdAbs. Assessing the Q-value on a per-residue level enabled us to identify residues that contribute to the instability and thus demonstrate which residues could be mutated to improve the stability and have later been validated by experiments. Our protocol extends beyond the application on sdAbs, as it is also suitable for other proteins and to determine the interfacial stability between protein and ligand.

Non-Ewald methods for evaluating the electrostatic interactions of charge systems: similarity and difference

In molecular simulations, it is essential to properly calculate the electrostatic interactions of particles in the physical system of interest. Here we consider a method called the non-Ewald method, which does not rely on the standard Ewald method with periodic boundary conditions, but instead relies on the cutoff-based techniques. We focus on the physicochemical and mathematical conceptual aspects of the method in order to gain a deeper understanding of the simulation methodology. In particular, we take into account the reaction field (RF) method, the isotropic periodic sum (IPS) method, and the zero-multipole summation method (ZMM). These cutoff-based methods are based on different physical ideas and are completely distinguishable in their underlying concepts. The RF and IPS methods are "additive" methods that incorporate information outside the cutoff region, via dielectric medium and isotropic boundary condition, respectively. In contrast, the ZMM is a "subtraction" method that tries to remove the artificial effects, generated near the boundary, from the cutoff sphere. Nonetheless, we find physical and/or mathematical similarities between these methods. In particular, the modified RF method can be derived by the principle of neutralization utilized in the ZMM, and we also found a direct relationship between IPS and ZMM.

Fly casting with ligand sliding and orientational selection supporting complex formation of a GPCR and a middle sized flexible molecule

A GA-guided multidimensional virtual-system coupled molecular dynamics (GA-mD-VcMD) simulation was conducted to elucidate binding mechanisms of a middle-sized flexible molecule, bosentan, to a GPCR protein, human endothelin receptor type B (hETB). GA-mD-VcMD is a generalized ensemble method that produces a free-energy landscape of the ligand-receptor binding by searching large-scale motions accompanied with stable maintenance of the fragile cell-membrane structure. All molecular components (bosentan, hETB, membrane, and solvent) were represented with an all-atom model. Then sampling was conducted from conformations where bosentan was distant from the binding site in the hETB binding pocket. The deepest basin in the resultant free-energy landscape was assigned to native-like complex conformation. The following binding mechanism was inferred. First, bosentan fluctuating randomly in solution is captured using a tip region of the flexible N-terminal tail of hETB via nonspecific attractive interactions (fly casting). Bosentan then slides occasionally from the tip to the root of the N-terminal tail (ligand–sliding). During this sliding, bosentan passes the gate of the binding pocket from outside to inside of the pocket with an accompanying rapid reduction of the molecular orientational variety of bosentan (orientational selection). Last, in the pocket, ligand–receptor attractive native contacts are formed. Eventually, the native-like complex is completed. The bosentan-captured conformations by the tip-region and root-region of the N-terminal tail correspond to two basins in the free-energy landscape. The ligand-sliding corresponds to overcoming of a free-energy barrier between the basins.

One-component plasma of a million particles via angular-averaged Ewald potential: A Monte Carlo study

In this work we derive a correct expression for the one-component plasma (OCP) energy via the angular-averaged Ewald potential (AAEP). Unlike Yakub and Ronchi [J. Low Temp. Phys. 139, 633 (2005)0022-229110.1007/s10909-005-5451-5], who had tried to obtain the same energy expression from a two-component plasma model, we used the original Ewald potential for an OCP. A constant in the AAEP was determined using the cluster expansion in the limit of weak coupling. The potential has a simple form suitable for effective numerical simulations. To demonstrate the advantages of the AAEP, we performed a number of Monte Carlo simulations for an OCP with up to a million particles in a wide range of the coupling parameter. Our computations turned out at least two orders of magnitude more effective than those with a traditional Ewald potential. A unified approach is offered for the determination of the thermodynamic limit in the whole investigated range. Our results are in good agreement with both theoretical data for a weakly coupled OCP and previous numerical simulations. We hope that the AAEP will be useful in path integral Monte Carlo simulations of the uniform electron gas.

Sphingolipids with 2-hydroxy fatty acids aid in plasma membrane nanodomain organization and oxidative burst

Plant sphingolipids mostly possess 2-hydroxy fatty acids (HFA), the synthesis of which is catalyzed by FA 2-hydroxylases (FAHs). In Arabidopsis (Arabidopsis thaliana), two FAHs (FAH1 and FAH2) have been identified. However, the functions of FAHs and sphingolipids with HFAs (2-hydroxy sphingolipids) are still unknown because of the lack of Arabidopsis lines with the complete deletion of FAH1. In this study, we generated a FAH1 mutant (fah1c) using CRISPR/Cas9-based genome editing. Sphingolipid analysis of fah1c, fah2, and fah1cfah2 mutants revealed that FAH1 hydroxylates very long-chain FAs (VLCFAs), whereas the substrates of FAH2 are VLCFAs and palmitic acid. However, 2-hydroxy sphingolipids are not completely lost in the fah1cfah2 double mutant, suggesting the existence of other enzymes catalyzing the hydroxylation of sphingolipid FAs. Plasma membrane (PM) analysis and molecular dynamics simulations revealed that hydroxyl groups of sphingolipid acyl chains play a crucial role in the organization of nanodomains, which are nanoscale liquid-ordered domains mainly formed by sphingolipids and sterols in the PM, through hydrogen bonds. In the PM of the fah1cfah2 mutant, the expression levels of 26.7% of the proteins, including defense-related proteins such as the pattern recognition receptors (PRRs) brassinosteroid insensitive 1-associated receptor kinase 1 and chitin elicitor receptor kinase 1, NADPH oxidase respiratory burst oxidase homolog D (RBOHD), and heterotrimeric G proteins, were lower than that in the wild-type. In addition, reactive oxygen species (ROS) burst was suppressed in the fah1cfah2 mutant after treatment with the pathogen-associated molecular patterns flg22 and chitin. These results indicated that 2-hydroxy sphingolipids are necessary for the organization of PM nanodomains and ROS burst through RBOHD and PRRs during pattern-triggered immunity.

N-Terminal-Driven Binding Mechanism of an Antigen Peptide to Human Leukocyte Antigen-A*2402 Elucidated by Multicanonical Molecular Dynamic-Based Dynamic Docking and Path Sampling Simulations

We have applied our advanced multicanonical molecular dynamics (McMD)-based dynamic docking methodology to investigate the binding mechanism of an HIV-1 Nef protein epitope to the Asian-dominant allele human leukocyte antigen (HLA)-A*2402. Even though pMHC complex formation [between a Major histocompatibility complex (MHC) class I molecule, which is encoded by an HLA allele, and an antigen peptide] is one of the fundamental processes of the adaptive human immune response, its binding mechanism has not yet been well studied, partially due to the high allelic variation of HLAs in the population. We have used our developed McMD-based dynamic docking method and have successfully reproduced the native complex structure, which is located near the free energy global minimum. Subsequent path sampling MD simulations elucidated the atomic details of the binding process and indicated that the peptide binding is initially driven by the highly positively charged N-terminus of the peptide that is attracted to the various negatively charged residues on the MHC molecule's surface. Upon nearing the pocket, the second tyrosine residue of the peptide anchors the peptide by strongly binding to the B-site of the MHC molecule via hydrophobic driven interactions, resulting in a very strong bound complex structure. Our methodology can be effectively used to predict the bound complex structures between MHC molecules and their antigens to study their binding mechanism in close detail, which would help with the development of new vaccines against cancers, as well as viral infections such as HIV and COVID-19.

Accurate Binding Configuration Prediction of a G-Protein-Coupled Receptor to Its Antagonist Using Multicanonical Molecular Dynamics-Based Dynamic Docking

We have performed dynamic docking between a prototypic G-protein-coupled receptor (GPCR) system, the β₂-adrenergic receptor, and its antagonist, alprenolol, using one of the enhanced conformation sampling methods, multicanonical molecular dynamics (McMD), which does not rely on any prior knowledge for the definition of the reaction coordinate. Although we have previously applied our McMD-based dynamic docking protocol to various globular protein systems, its application to GPCR systems would be difficult because of their complicated design, which include a lipid bilayer, and because of the difficulty in sampling the configurational space of a binding site that exists deep inside the GPCR. Our simulations sampled a wide array of ligand-bound and ligand-unbound structures, and we measured 427 binding events during our 48 μs production run. Analysis of the ensemble revealed several stable and meta-stable structures, where the most stable structure at the global free energy minimum matches the experimental one. Additional canonical MD simulations were used for refinement and validation of the structures, revealing that most of the intermediates are sufficiently stable to trap the ligand in these intermediary states and furthermore validated our prediction results. Given the difficulty in reaching the orthosteric binding site, chemical optimization of the compound for the second ranking configuration, which binds near the pocket's entrance, might lead to a high-affinity allosteric inhibitor. Accordingly, we show that the application of our methodology can be used to provide crucial insights for the rational design of drugs that target GPCRs.

Enzyme modification using mutation site prediction method for enhancing the regioselectivity of substrate reaction sites

Enzymes with low regioselectivity of substrate reaction sites may produce multiple products from a single substrate. When a target product is produced industrially using these enzymes, the production of non-target products (byproducts) causes adverse effects such as increased processing costs for purification and the amount of raw material. Thus it is required the development of modified enzymes to reduce the amount of byproducts’ production. In this paper, we report a method called mutation site prediction for enhancing the regioselectivity of substrate reaction sites (MSPER). MSPER takes conformational data for docking poses of an enzyme and a substrate as input and automatically generates a ranked list of mutation sites to destabilize docking poses for byproducts while maintaining those for target products in silico. We applied MSPER to the enzyme cytochrome P450 CYP102A1 (BM3) and the two substrates to enhance the regioselectivity for four target products with different reaction sites. The 13 of the total 14 top-ranked mutation sites predicted by MSPER for the four target products succeeded in selectively enhancing the regioselectivity up to 6.4-fold. The results indicate that MSPER can distinguish differences of substrate structures and the reaction sites, and can accurately predict mutation sites to enhance regioselectivity without selection by directed evolution screening.

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.

Difference of binding modes among three ligands to a receptor mSin3B corresponding to their inhibitory activities

A preceding experiment suggested that a compound, which inhibits binding of the REST/NRSF segment to the cleft of a receptor protein mSin3B, can be a potential drug candidate to ameliorate many neuropathies. We have recently developed an enhanced conformational sampling method, genetic-algorithm-guided multi-dimensional virtual-system-coupled canonical molecular dynamics, and in the present study, applied it to three systems consisting of mSin3B and one of three compounds, sertraline, YN3, and acitretin. Other preceding experiments showed that only sertraline inhibits the binding of REST/NRSF to mSin3B. The current simulation study produced the spatial distribution of the compounds around mSin3B, and showed that sertraline and YN3 bound to the cleft of mSin3B with a high propensity, although acitretin did not. Further analyses of the simulation data indicated that only the sertraline-mSin3B complex produced a hydrophobic core similar to that observed in the molecular interface of the REST/NRSF-mSin3B complex: An aromatic ring of sertraline sunk deeply in the mSin3B's cleft forming a hydrophobic core contacting to hydrophobic amino-acid residues located at the bottom of the cleft. The present study proposes a step to design a compound that inhibits competitively the binding of a ligand to its receptor.

Cryptic-site binding mechanism of medium-sized Bcl-xL inhibiting compounds elucidated by McMD-based dynamic docking simulations

We have performed multicanonical molecular dynamics (McMD) based dynamic docking simulations to study and compare the binding mechanism between two medium-sized inhibitors (ABT-737 and WEHI-539) that bind to the cryptic site of Bcl-xL, by exhaustively sampling the conformational and configurational space. Cryptic sites are binding pockets that are transiently formed in the apo state or are induced upon ligand binding. Bcl-xL, a pro-survival protein involved in cancer progression, is known to have a cryptic site, whereby the shape of the pocket depends on which ligand is bound to it. Starting from the apo-structure, we have performed two independent McMD-based dynamic docking simulations for each ligand, and were able to obtain near-native complex structures in both cases. In addition, we have also studied their interactions along their respective binding pathways by using path sampling simulations, which showed that the ligands form stable binding configurations via predominantly hydrophobic interactions. Although the protein started from the apo state, both ligands modulated the pocket in different ways, shifting the conformational preference of the sub-pockets of Bcl-xL. We demonstrate that McMD-based dynamic docking is a powerful tool that can be effectively used to study binding mechanisms involving a cryptic site, where ligand binding requires a large conformational change in the protein to occur.

Dielectric continuum model examination of real-space electrostatic treatments

Electrostatic interaction is long ranged; thus, the accurate calculation is not an easy task in molecular dynamics or Monte Carlo simulations. Though the rigorous Ewald method based on the reciprocal space has been established, real-space treatments have recently become an attractive alternative because of the efficient calculation. However, the construction is not yet completed and is now a challenging subject. In an earlier theoretical study, Neumann and Steinhauser employed the Onsager dielectric continuum model to explain how simple real-space cutoff produces artificial dipolar orientation. In the present study, we employ this continuum model to explore the fundamental properties of the recently developed real-space treatments of three shifting schemes. The result of the distance-dependent Kirkwood function G_K(R) showed that the simple bare cutoff produces a well-known hole-shaped artifact, whereas the shift treatments do not. Two-dimensional mapping of electric field well explained how these shift treatments remove the hole-shaped artifact. Still, the shift treatments are not sufficient because they do not produce a flat G_K(R) profile unlike ideal no-cutoff treatment. To test the continuum model results, we also performed Monte Carlo simulations of dipolar particles. The results found that the continuum model could predict the qualitative tendency as to whether each electrostatic treatment produces the hole-shaped artifact of G_K(R) or not. We expect that the present study using the continuum model offers a stringent criterion to judge whether the primitive electrostatic behavior is correctly described or not, which will be useful for future construction of electrostatic treatments.

Dynamic Docking Using Multicanonical Molecular Dynamics: Simulating Complex Formation at the Atomistic Level

Multicanonical molecular dynamics (McMD)-based dynamic docking has been applied to predict the native binding configurations for several protein receptors and their ligands. Due to the enhanced sampling capabilities of McMD, it can exhaustively sample bound and unbound ligand configurations, as well as receptor conformations, and thus enables efficient sampling of the conformational and configurational space, not possible using canonical MD simulations. As McMD samples a wide configurational space, extensive analysis is required to study the diverse ensemble consisting of bound and unbound structures. By projecting the reweighted ensemble onto the first two principal axes obtained via principal component analysis of the multicanonical ensemble, the free energy landscape (FEL) can be obtained. Further analysis produces representative structures positioned at the local minima of the FEL, where these structures are then ranked by their free energy. In this chapter, we describe our dynamic docking methodology, which has successfully reproduced the native binding configuration for small compounds, medium-sized compounds, and peptide molecules.

Assessment of the Wolf method using the Stillinger-Lovett sum rules: From strong electrolytes to weakly charged colloidal dispersions

The Ewald method has been the cornerstone in molecular simulations for modeling electrostatic interactions of charge-stabilized many-body systems. In the late 1990s, Wolf and collaborators developed an alternative route to describe the long-range nature of electrostatic interactions; from a computational perspective, this method provides a more efficient and straightforward way to implement long-range electrostatic interactions than the Ewald method. Despite these advantages, the validity of the Wolf potential to account for the electrostatic contribution in charged fluids remains controversial. To alleviate this situation, in this contribution, we implement the Wolf summation method to both electrolyte solutions and charged colloids with moderate size and charge asymmetries in order to assess the accuracy and validity of the method. To this end, we verify that the proper selection of parameters within the Wolf method leads to results that are in good agreement with those obtained through the standard Ewald method and the theory of integral equations of simple liquids within the so-called hypernetted chain approximation. Furthermore, we show that the results obtained with the original Wolf method do satisfy the moment conditions described by the Stillinger-Lovett sum rules, which are directly related to the local electroneutrality condition and the electrostatic screening in the Debye-Hückel regime. Hence, the fact that the solution provided by the Wolf method satisfies the first and second moments of Stillinger-Lovett proves, for the first time, the reliability of the method to correctly incorporate the electrostatic contribution in charge-stabilized fluids. This makes the Wolf method a powerful alternative compared to more demanding computational approaches.

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ε.

Exhaustive search of the configurational space of heat-shock protein 90 with its inhibitor by multicanonical molecular dynamics based dynamic docking

Multicanonical molecular dynamics based dynamic docking was used to exhaustively search the configurational space of an inhibitor binding to the N-terminal domain of heat-shock protein 90 (Hsp90). The obtained structures at 300 K cover a wide structural ensemble, with the top two clusters ranked by their free energy coinciding with the native binding site. The representative structure of the most stable cluster reproduced the experimental binding configuration, but an interesting conformational change in Hsp90 could be observed. The combined effects of solvation and ligand binding shift the equilibrium from a preferred loop-in conformation in the unbound state to an α-helical one in the bound state for the flexible lid region of Hsp90. Thus, our dynamic docking method is effective at predicting the native binding site while exhaustively sampling a wide configurational space, modulating the protein structure upon binding.

Generalized Moment Correction for Long-Ranged Electrostatics

Describing long-ranged electrostatics using short-ranged pair potentials is appealing because the computational complexity scales linearly with the number of particles. The foundation of the approach presented here is to mimic the long-ranged medium response by cancelling electric multipoles within a small cutoff sphere. We propose a rigorous and formally exact new method that cancels up to infinitely many multipole moments and is free of operational damping parameters often required in existing theories. Using molecular dynamics simulations of water with and without added salt, we discuss radial distribution functions, Kirkwood–Buff integrals, dielectrics, diffusion coefficients, and angular correlations in relation to existing electrostatic models. We find that the proposed method is an efficient and accurate alternative for handling long-ranged electrostatics as compared to Ewald summation schemes. The methodology and proposed parameterization are applicable also for dipole–dipole interactions.

Reaction-field electrostatics in molecular dynamics simulations: development of a conservative scheme compatible with an atomic cutoff

Shifting and switching schemes are developed to enable strict energy conservation in molecular dynamics simulations relying on reaction-field electrostatic (as well as Lennard-Jones) interactions with an atom-based cutoff truncation.

Free-energy landscape of molecular interactions between endothelin 1 and human endothelin type B receptor: fly-casting mechanism

The free-energy landscape of interaction between a medium-sized peptide, endothelin 1 (ET1), and its receptor, human endothelin type B receptor (hETB), was computed using multidimensional virtual-system coupled molecular dynamics, which controls the system's motions by introducing multiple reaction coordinates. The hETB embedded in lipid bilayer was immersed in explicit solvent. All molecules were expressed as all-atom models. The resultant free-energy landscape had five ranges with decreasing ET1-hETB distance: completely dissociative, outside-gate, gate, binding pocket, and genuine-bound ranges. In the completely dissociative range, no ET1-hETB interaction appeared. In the outside-gate range, an ET1-hETB attractive interaction was the fly-casting mechanism. In the gate range, the ET1 orientational variety decreased rapidly. In the binding pocket range, ET1 was in a narrow pathway with a steep free-energy slope. In the genuine-bound range, ET1 was in a stable free-energy basin. A G-protein-coupled receptor (GPCR) might capture its ligand from a distant place.

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.

Effects of number of parallel runs and frequency of bias-strength replacement in generalized ensemble molecular dynamics simulations

Background

The generalized ensemble approach with the molecular dynamics (MD) method has been widely utilized. This approach usually has two features. (i) A bias potential, whose strength is replaced during a simulation, is applied. (ii) Sampling can be performed by many parallel runs of simulations. Although the frequency of the bias-strength replacement and the number of parallel runs can be adjusted, the effects of these settings on the resultant ensemble remain unclear.

Method

In this study, we performed multicanonical MD simulations for a foldable mini-protein (Trp-cage) and two unstructured peptides (8- and 20-residue poly-glutamic acids) with various settings.

Results

As a result, running many short simulations yielded robust results for the Trp-cage model. Regarding the frequency of the bias-potential replacement, although using a high frequency enhanced the traversals in the potential energy space, it did not promote conformational changes in all the systems.

Theoretical equations of state for a charged fluid

In this article, we present a molecular thermodynamic study of a system of N particles contained within a volume V and interacting via a hard-core pair potential with an attractive interaction according to the Wolf model for charged systems. This variable-range potential is characterized by three parameters: the repulsive hard-core diameter σ, the energy-well depth ϵ, and the inverse range α; a fourth parameter of the model is a cut-off distance x_c that depends on α according to the relation x_c = 2/α. Two equations of state (EOSs) are presented and derived from thermodynamic perturbation theory and Monte Carlo (MC) simulation data. The first EOS is given by the standard Zwanzig's high-temperature expansion of the Helmholtz free energy, where the first three perturbation terms a₁, a₂, and a₃ were obtained from MC simulations in the canonical ensemble (NVT) and parameterized as functions of α and the reduced density of particles ρ^* = Nσ³/V. The second EOS was obtained from the discrete perturbation theory applied to a discrete representation of the Wolf potential. Results for pressures, internal energies, and isochoric heat capacities are compared to the MC computer simulation data of the Wolf system, including vapor-liquid coexistence curves, for different values of α. Overall, both EOSs give a very good representation of the thermodynamic properties of the Wolf fluid when 0.3 ≤ α ≤ 1.0 and 0.05 ≤ ρ^* ≤ 0.8. Since the Yukawa fluid can reproduce information of screened ionic interactions, we discuss the equivalence between the Wolf and Yukawa fluids in the context of equivalent systems in liquid theory.

On short-ranged pair-potentials for long-range electrostatics

Fast and accurate summation of long-range electrostatics by using a short-ranged pair-potential that ensures moment cancellation in the cutoff sphere.

A Modified Shifted Force Approach to the Wolf Summation

The Wolf method for calculation of electrostatic interactions in molecular simulations is known to describe the energy well, whereas the forces have discontinuities. For a more reliable description of the forces this method can be extended with a shifted force approach. This leads to a good description of the forces and precise molecular dynamics simulation, but the description of the energy becomes poorer. In this study we propose a modification of a shifted force extension to describe the energy as well as the forces in better agreement to reference data as determined from the Ewald summation. We show that vapor-liquid phase equilibria (VLE) calculated with Monte Carlo simulations in the grand canonical ensemble and dynamic properties calculated with molecular dynamics simulations can be calculated reliably using this modification to describe the electrostatic interactions.

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.

Enhanced Sampling of Molecular Dynamics Simulations of a Polyalanine Octapeptide: Effects of the Periodic Boundary Conditions on Peptide Conformation

We investigate the problem of artifacts caused by the periodic boundary conditions (PBC) used in molecular simulation studies. Despite the long history of simulations with PBCs, the existence of measurable artifacts originating from PBCs applied to inherently nonperiodic physical systems remains controversial. Specifically, these artifacts appear as differences between simulations of the same system but with different simulation-cell sizes. Earlier studies have implied that, even in the simple case of a small model peptide in water, sampling inefficiency is a major obstacle to understanding these artifacts. In this study, we have resolved the sampling issue using the replica exchange molecular dynamics (REMD) enhanced-sampling method to explore PBC artifacts. Explicitly solvated zwitterionic polyalanine octapeptides with three different cubic-cells, having dimensions of L = 30, 40, and 50 Å, were investigated to elucidate the differences with 64 replica × 500 ns REMD simulations using the AMBER parm99SB force field. The differences among them were not large overall, and the results for the L = 30 and 40 Å simulations in the conformational free energy landscape were found to be very similar at room temperature. However, a small but statistically significant difference was seen for L = 50 Å. We observed that extended conformations were slightly overstabilized in the smaller systems. The origin of these artifacts is discussed by comparison to an electrostatic calculation method without PBCs.

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.

Molecular mechanisms of cooperative binding of transcription factors Runx1-CBFβ-Ets1 on the TCRα gene enhancer

Ets1 is an essential transcription factor (TF) for several important physiological processes, including cell proliferation and differentiation. Its recognition of the enhancer region of the TCRα gene is enhanced by the cooperative binding of the Runx1–CBFβ heterodimer, with the cancelation of phosphorylation-dependent autoinhibition. The detailed mechanism of this interesting cooperativity between Ets1 and the Runx1–CBFβ heterodimer is still largely unclear. Here, we investigated the molecular mechanisms of this cooperativity, by using molecular dynamics simulations. Consequently, we detected high flexibility of the loop region between the HI2 and H1 helices of Ets1. Upon Runx1–CBFβ heterodimer binding, this loop transiently adopts various sub-stable conformations in its interactions with the DNA. In addition, a network analysis suggested an allosteric pathway in the molecular assembly and identified some key residues that coincide with previous experimental studies. Our simulations suggest that the cooperative binding of Ets1 and the Runx1–CBFβ heterodimer alters the DNA conformation and induces sub-stable conformations of the HI2–H1 loop of Ets1. This phenomenon increases the flexibility of the regulatory module, including the HI2 helix, and destabilizes the inhibitory form of this module. Thus, we hypothesize that this effect facilitates Ets1–DNA binding and prevents the phosphorylation-dependent DNA binding autoinhibition.

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.

Revisiting antibody modeling assessment for CDR-H3 loop

The antigen-binding site of antibodies, also known as complementarity-determining region (CDR), has hypervariable sequence properties. In particular, the third CDR loop of the heavy chain, CDR-H3, has such variability in its sequence, length, and conformation that ordinary modeling techniques cannot build a high-quality structure. At Stage 2 of the Second Antibody Modeling Assessment (AMA-II) held in 2013, the model structures of the CDR-H3 loops were submitted by the seven modelers and were critically assessed. After our participation in AMA-II, we rebuilt one of the long CDR-H3 loops with 13 residues (A52 antibody) by a more precise method, using enhanced conformational sampling with the explicit water model, as compared to our previous method employed at AMA-II. The current stable models obtained from the free energy landscape at 300 K include structures similar to the X-ray crystal structures. Those models were not built in our previous work at AMA-II. The current free energy landscape suggested that the CDR-H3 loop structures in the crystal are not stable in solution, but they are stabilized by the crystal packing effect.

myPresto/omegagene: a GPU-accelerated molecular dynamics simulator tailored for enhanced conformational sampling methods with a non-Ewald electrostatic scheme

Molecular dynamics (MD) is a promising computational approach to investigate dynamical behavior of molecular systems at the atomic level. Here, we present a new MD simulation engine named "myPresto/omegagene" that is tailored for enhanced conformational sampling methods with a non-Ewald electrostatic potential scheme. Our enhanced conformational sampling methods, e.g., the virtual-system-coupled multi-canonical MD (V-McMD) method, replace a multi-process parallelized run with multiple independent runs to avoid inter-node communication overhead. In addition, adopting the non-Ewald-based zero-multipole summation method (ZMM) makes it possible to eliminate the Fourier space calculations altogether. The combination of these state-of-the-art techniques realizes efficient and accurate calculations of the conformational ensemble at an equilibrium state. By taking these advantages, myPresto/omegagene is specialized for the single process execution with Graphics Processing Unit (GPU). We performed benchmark simulations for the 20-mer peptide, Trp-cage, with explicit solvent. One of the most thermodynamically stable conformations generated by the V-McMD simulation is very similar to an experimentally solved native conformation. Furthermore, the computation speed is four-times faster than that of our previous simulation engine, myPresto/psygene-G. The new simulator, myPresto/omegagene, is freely available at the following URLs: http://www.protein.osaka-u.ac.jp/rcsfp/pi/omegagene/ and http://presto.protein.osaka-u.ac.jp/myPresto4/.

Elastic properties of dynein motor domain obtained from all-atom molecular dynamics simulations

Dyneins are large microtubule motor proteins that convert ATP energy to mechanical power. High-resolution crystal structures of ADP-bound cytoplasmic dynein have revealed the organization of the motor domain, comprising the AAA(+) ring, the linker, the stalk/strut and the C sequence. Recently, the ADP.vanadate-bound structure, which is similar to the ATP hydrolysis transition state, revealed how the structure of dynein changes upon ATP binding. Although both the ADP- and ATP-bound state structures have been resolved, the dynamic properties at the atomic level remain unclear. In this work, we built two models named 'the ADP model' and 'the ATP model', where ADP and ATP are bound to AAA1 in the AAA(+) ring, respectively, to observe the initial procedure of the structural change from the unprimed to the primed state. We performed 200-ns molecular dynamics simulations for both models and compared their structures and dynamics. The motions of the stalk, consisting of a long coiled coil with a microtubule-binding domain, significantly differed between the two models. The elastic properties of the stalk were analyzed and compared with the experimental results.