A revised density function for molecular surface definition in continuum solvent models

Calculation of solvation force in molecular dynamics simulation by deep-learning method

Electrostatic calculations are generally used in studying the thermodynamics and kinetics of biomolecules in solvent. Generally, this is performed by solving the Poisson-Boltzmann equation on a large grid system, a process known to be time consuming. In this study, we developed a deep neural network to predict the decomposed solvation free energies and forces of all atoms in a molecule. To train the network, the internal coordinates of the molecule were used as the input data, and the solvation free energies along with transformed atomic forces from the Poisson-Boltzmann equation were used as labels. Both the training and prediction tasks were accelerated on GPU. Formal tests demonstrated that our method can provide reasonable predictions for small molecules when the network is well-trained with its simulation data. This method is suitable for processing lots of snapshots of molecules in a long trajectory. Moreover, we applied this method in the molecular dynamics simulation with enhanced sampling. The calculated free energy landscape closely resembled that obtained from explicit solvent simulations.

Grid-Robust Efficient Neural Interface Model for Universal Molecule Surface Construction from Point Clouds

Molecular surfaces play a pivotal role in elucidating the properties and functions of biological complexes. While various surfaces have been proposed for specific scenarios, their widespread adoption faces challenges due to limited efficiency stemming from hand-crafted modeling designs. In this work, we proposed a general framework that incorporates both the point cloud concept and neural networks. The use of matrix multiplication in this framework enables efficient implementation across diverse platforms and libraries. We applied this framework to develop the GENIUSES (Grid-robust Efficient Neural Interface for Universal Solvent-Excluded Surface) model for constructing SES. GENIUSES demonstrates high accuracy and efficiency across data sets with varying conformations and complexities. Compared to the classical implementation of SES in the AMBER software package, our framework achieved a 26-fold speedup while retaining ∼95% accuracy when ported to the GPU platform using CUDA. Greater speedups can be obtained in large-scale systems. Importantly, our model exhibits robustness against variations in the grid spacing. We have integrated this infrastructure into AMBER to enhance accessibility for research in drug screening and related fields, where efficiency is of paramount importance.

Molecular dynamics study reveals key disruptors of MEIG1-PACRG interaction

Interactions between the meiosis-expressed gene 1 (MEIG1) and Parkin co-regulated gene (PACRG) protein are critical in the formation of mature sperm cells. Targeting either MEIG1 or PACRG protein could be a contraceptive strategy. The W50A and Y68A mutations on MEIG1 are known to interrupt the MEIG1-PACRG interactions resulting in defective sperm cells. However, the details about how the mutants disrupt the protein-protein binding are not clear. In this study, we reveal insights on MEIG1 and PACRG protein dynamics by applying Gaussian-accelerated molecular dynamics (GaMD) simulations and post-GaMD analysis. Our results show that the mutations destabilize the protein-protein interfacial interaction. The effect of the Y68A mutation is more significant than W50A as Y68 forms stronger polar interactions with PACRG. Because both human and mouse models demonstrate similar dynamic properties, the findings from mouse proteins can be applied to the human system. Moreover, we report a potential ligand binding pocket on the MEIG1 and PACRG interaction surface that could be a target for future drug design to inhibit the MEIG1-PACRG interaction. PACRG shows more qualified pockets along the protein-protein interface, implying that it is a better target than MEIG1. Our work provides a fundamental understanding of MEIG1 and PACRG protein dynamics, paving the way for drug discovery in male-based contraception.

Exploring the druggability of the binding site of aurovertin, an exogenous allosteric inhibitor of FOF1-ATP synthase

In addition to playing a central role in the mitochondria as the main producer of ATP, F_OF₁-ATP synthase performs diverse key regulatory functions in the cell membrane. Its malfunction has been linked to a growing number of human diseases, including hypertension, atherosclerosis, cancer, and some neurodegenerative, autoimmune, and aging diseases. Furthermore, inhibition of this enzyme jeopardizes the survival of several bacterial pathogens of public health concern. Therefore, F_OF₁-ATP synthase has emerged as a novel drug target both to treat human diseases and to combat antibiotic resistance. In this work, we carried out a computational characterization of the binding sites of the fungal antibiotic aurovertin in the bovine F₁ subcomplex, which shares a large identity with the human enzyme. Molecular dynamics simulations showed that although the binding sites can be described as preformed, the inhibitor hinders inter-subunit communications and exerts long-range effects on the dynamics of the catalytic site residues. End-point binding free energy calculations revealed hot spot residues for aurovertin recognition. These residues were also relevant to stabilize solvent sites determined from mixed-solvent molecular dynamics, which mimic the interaction between aurovertin and the enzyme, and could be used as pharmacophore constraints in virtual screening campaigns. To explore the possibility of finding species-specific inhibitors targeting the aurovertin binding site, we performed free energy calculations for two bacterial enzymes with experimentally solved 3D structures. Finally, an analysis of bacterial sequences was carried out to determine conservation of the aurovertin binding site. Taken together, our results constitute a first step in paving the way for structure-based development of new allosteric drugs targeting F_OF₁-ATP synthase sites of exogenous inhibitors.

Machine-Learned Molecular Surface and Its Application to Implicit Solvent Simulations

Implicit solvent models, such as Poisson-Boltzmann models, play important roles in computational studies of biomolecules. A vital step in almost all implicit solvent models is to determine the solvent-solute interface, and the solvent excluded surface (SES) is the most widely used interface definition in these models. However, classical algorithms used for computing SES are geometry-based, so that they are neither suitable for parallel implementations nor convenient for obtaining surface derivatives. To address the limitations, we explored a machine learning strategy to obtain a level set formulation for the SES. The training process was conducted in three steps, eventually leading to a model with over 95% agreement with the classical SES. Visualization of tested molecular surfaces shows that the machine-learned SES overlaps with the classical SES in almost all situations. Further analyses show that the machine-learned SES is incredibly stable in terms of rotational variation of tested molecules. Our timing analysis shows that the machine-learned SES is roughly 2.5 times as efficient as the classical SES routine implemented in Amber/PBSA on a tested central processing unit (CPU) platform. We expect further performance gain on massively parallel platforms such as graphics processing units (GPUs) given the ease in converting the machine-learned SES to a parallel procedure. We also implemented the machine-learned SES into the Amber/PBSA program to study its performance on reaction field energy calculation. The analysis shows that the two sets of reaction field energies are highly consistent with a 1% deviation on average. Given its level set formulation, we expect the machine-learned SES to be applied in molecular simulations that require either surface derivatives or high efficiency on parallel computing platforms.

Recent Developments in Free Energy Calculations for Drug Discovery

The grand challenge in structure-based drug design is achieving accurate prediction of binding free energies. Molecular dynamics (MD) simulations enable modeling of conformational changes critical to the binding process, leading to calculation of thermodynamic quantities involved in estimation of binding affinities. With recent advancements in computing capability and predictive accuracy, MD based virtual screening has progressed from the domain of theoretical attempts to real application in drug development. Approaches including the Molecular Mechanics Poisson Boltzmann Surface Area (MM-PBSA), Linear Interaction Energy (LIE), and alchemical methods have been broadly applied to model molecular recognition for drug discovery and lead optimization. Here we review the varied methodology of these approaches, developments enhancing simulation efficiency and reliability, remaining challenges hindering predictive performance, and applications to problems in the fields of medicine and biochemistry.

Improved Poisson-Boltzmann Methods for High-Performance Computing

Implicit solvent models based on the Poisson-Boltzmann equation (PBE) have been widely used to study electrostatic interactions in biophysical processes. These models often treat the solvent and solute regions as high and low dielectric continua, leading to a large jump in dielectrics across the molecular surface which is difficult to handle. Higher order interface schemes are often needed to seek higher accuracy for PBE applications. However, these methods are usually very liberal in the use of grid points nearby the molecular surface, making them difficult to use on high-performance computing platforms. Alternatively, the harmonic average (HA) method has been used to approximate dielectric interface conditions near the molecular surface with surprisingly good convergence and is well suited for high-performance computing. By adopting a 7-point stencil, the HA method is advantageous in generating simple 7-banded coefficient matrices, which greatly facilitate linear system solution with dense data parallelism, on high-performance computing platforms such as a graphics processing unit (GPU). However, the HA method is limited due to its lower accuracy. Therefore, it would be of great interest for high-performance applications to develop more accurate methods while retaining the simplicity and effectiveness of the 7-point stencil discretization scheme. In this study, we have developed two new algorithms based on the spirit of the HA method by introducing more physical interface relations and imposing the discretized Poisson's equation to the second order, respectively. Our testing shows that, for typical biomolecules, the new methods significantly improve the numerical accuracy to that comparable to the second-order solvers and with ∼65% overall efficiency gain on widely available high-performance GPU platforms.

Heterogeneous Dielectric Implicit Membrane Model for the Calculation of MMPBSA Binding Free Energies

Membrane-bound protein receptors are a primary biological drug target, but the computational analysis of membrane proteins has been limited. In order to improve molecular mechanics Poisson-Boltzmann surface area (MMPBSA) binding free energy calculations for membrane protein-ligand systems, we have optimized a new heterogeneous dielectric implicit membrane model, with respect to free energy simulations in explicit membrane and explicit water, and implemented it into the Amber software suite. This new model supersedes our previous uniform, single dielectric implicit membrane model by allowing the dielectric constant to vary with depth within the membrane. We calculated MMPBSA binding free energies for the human purinergic platelet receptor (P2Y₁₂R) and two of the muscarinic acetylcholine receptors (M2R and M3R) bound to various antagonist ligands using both membrane models, and we found that the heterogeneous dielectric membrane model has a stronger correlation with experimental binding affinities compared to the older model under otherwise identical conditions. This improved membrane model increases the utility of MMPBSA calculations for the rational design and improvement of future drug candidates.

An efficient second-order poisson-boltzmann method

Immersed interface method (IIM) is a promising high-accuracy numerical scheme for the Poisson-Boltzmann model that has been widely used to study electrostatic interactions in biomolecules. However, the IIM suffers from instability and slow convergence for typical applications. In this study, we introduced both analytical interface and surface regulation into IIM to address these issues. The analytical interface setup leads to better accuracy and its convergence closely follows a quadratic manner as predicted by theory. The surface regulation further speeds up the convergence for nontrivial biomolecules. In addition, uncertainties of the numerical energies for tested systems are also reduced by about half. More interestingly, the analytical setup significantly improves the linear solver efficiency and stability by generating more precise and better-conditioned linear systems. Finally, we implemented the bottleneck linear system solver on GPUs to further improve the efficiency of the method, so it can be widely used for practical biomolecular applications. © 2019 Wiley Periodicals, Inc.

Robustness and Efficiency of Poisson-Boltzmann Modeling on Graphics Processing Units

Poisson-Boltzmann equation (PBE) based continuum electrostatics models have been widely used in modeling electrostatic interactions in biochemical processes, particularly in estimating protein-ligand binding affinities. Fast convergence of PBE solvers is crucial in binding affinity computations as numerous snapshots need to be processed. Efforts have been reported to develop PBE solvers on graphics processing units (GPUs) for efficient modeling of biomolecules, though only relatively simple successive over-relaxation and conjugate gradient methods were implemented. However, neither convergence nor scaling properties of the two methods are optimal for large biomolecules. On the other hand, geometric multigrid (MG) has been shown to be an optimal solver on CPUs, though no MG have been reported for biomolecular applications on GPUs. This is not a surprise as it is a more complex method and depends on simpler but limited iterative methods such as Gauss-Seidel in its core relaxation procedure. The robustness and efficiency of MG on GPUs are also unclear. Here we present an implementation and a thorough analysis of MG on GPUs. Our analysis shows that robustness is a more pronounced issue than efficiency for both MG and other tested solvers when the single precision is used for complex biomolecules. We further show how to balance robustness and efficiency utilizing MG's overall efficiency and conjugate gradient's robustness, pointing to a hybrid GPU solver with a good balance of efficiency and accuracy. The new PBE solver will significantly improve the computational throughput for a range of biomolecular applications on the GPU platforms.

Exploring a multi-scale method for molecular simulation in continuum solvent model: Explicit simulation of continuum solvent as an incompressible fluid

We explored a multi-scale algorithm for the Poisson-Boltzmann continuum solvent model for more robust simulations of biomolecules. In this method, the continuum solvent/solute interface is explicitly simulated with a numerical fluid dynamics procedure, which is tightly coupled to the solute molecular dynamics simulation. There are multiple benefits to adopt such a strategy as presented below. At this stage of the development, only nonelectrostatic interactions, i.e., van der Waals and hydrophobic interactions, are included in the algorithm to assess the quality of the solvent-solute interface generated by the new method. Nevertheless, numerical challenges exist in accurately interpolating the highly nonlinear van der Waals term when solving the finite-difference fluid dynamics equations. We were able to bypass the challenge rigorously by merging the van der Waals potential and pressure together when solving the fluid dynamics equations and by considering its contribution in the free-boundary condition analytically. The multi-scale simulation method was first validated by reproducing the solute-solvent interface of a single atom with analytical solution. Next, we performed the relaxation simulation of a restrained symmetrical monomer and observed a symmetrical solvent interface at equilibrium with detailed surface features resembling those found on the solvent excluded surface. Four typical small molecular complexes were then tested, both volume and force balancing analyses showing that these simple complexes can reach equilibrium within the simulation time window. Finally, we studied the quality of the multi-scale solute-solvent interfaces for the four tested dimer complexes and found that they agree well with the boundaries as sampled in the explicit water simulations.

Recent Developments and Applications of the MMPBSA Method

The Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) approach has been widely applied as an efficient and reliable free energy simulation method to model molecular recognition, such as for protein-ligand binding interactions. In this review, we focus on recent developments and applications of the MMPBSA method. The methodology review covers solvation terms, the entropy term, extensions to membrane proteins and high-speed screening, and new automation toolkits. Recent applications in various important biomedical and chemical fields are also reviewed. We conclude with a few future directions aimed at making MMPBSA a more robust and efficient method.

A Continuum Poisson-Boltzmann Model for Membrane Channel Proteins

Membrane proteins constitute a large portion of the human proteome and perform a variety of important functions as membrane receptors, transport proteins, enzymes, signaling proteins, and more. Computational studies of membrane proteins are usually much more complicated than those of globular proteins. Here, we propose a new continuum model for Poisson-Boltzmann calculations of membrane channel proteins. Major improvements over the existing continuum slab model are as follows: (1) The location and thickness of the slab model are fine-tuned based on explicit-solvent MD simulations. (2) The highly different accessibilities in the membrane and water regions are addressed with a two-step, two-probe grid-labeling procedure. (3) The water pores/channels are automatically identified. The new continuum membrane model is optimized (by adjusting the membrane probe, as well as the slab thickness and center) to best reproduce the distributions of buried water molecules in the membrane region as sampled in explicit water simulations. Our optimization also shows that the widely adopted water probe of 1.4 Å for globular proteins is a very reasonable default value for membrane protein simulations. It gives the best compromise in reproducing the explicit water distributions in membrane channel proteins, at least in the water accessible pore/channel regions. Finally, we validate the new membrane model by carrying out binding affinity calculations for a potassium channel, and we observe good agreement with the experimental results.

Acceleration of Linear Finite-Difference Poisson-Boltzmann Methods on Graphics Processing Units

Electrostatic interactions play crucial roles in biophysical processes such as protein folding and molecular recognition. Poisson-Boltzmann equation (PBE)-based models have emerged as widely used in modeling these important processes. Though great efforts have been put into developing efficient PBE numerical models, challenges still remain due to the high dimensionality of typical biomolecular systems. In this study, we implemented and analyzed commonly used linear PBE solvers for the ever-improving graphics processing units (GPU) for biomolecular simulations, including both standard and preconditioned conjugate gradient (CG) solvers with several alternative preconditioners. Our implementation utilizes the standard Nvidia CUDA libraries cuSPARSE, cuBLAS, and CUSP. Extensive tests show that good numerical accuracy can be achieved given that the single precision is often used for numerical applications on GPU platforms. The optimal GPU performance was observed with the Jacobi-preconditioned CG solver, with a significant speedup over standard CG solver on CPU in our diversified test cases. Our analysis further shows that different matrix storage formats also considerably affect the efficiency of different linear PBE solvers on GPU, with the diagonal format best suited for our standard finite-difference linear systems. Further efficiency may be possible with matrix-free operations and integrated grid stencil setup specifically tailored for the banded matrices in PBE-specific linear systems.

Numerical interpretation of molecular surface field in dielectric modeling of solvation

Continuum solvent models, particularly those based on the Poisson-Boltzmann equation (PBE), are widely used in the studies of biomolecular structures and functions. Existing PBE developments have been mainly focused on how to obtain more accurate and/or more efficient numerical potentials and energies. However to adopt the PBE models for molecular dynamics simulations, a difficulty is how to interpret dielectric boundary forces accurately and efficiently for robust dynamics simulations. This study documents the implementation and analysis of a range of standard fitting schemes, including both one-sided and two-sided methods with both first-order and second-order Taylor expansions, to calculate molecular surface electric fields to facilitate the numerical calculation of dielectric boundary forces. These efforts prompted us to develop an efficient approximated one-dimensional method, which is to fit the surface field one dimension at a time, for biomolecular applications without much compromise in accuracy. We also developed a surface-to-atom force partition scheme given a level set representation of analytical molecular surfaces to facilitate their applications to molecular simulations. Testing of these fitting methods in the dielectric boundary force calculations shows that the second-order methods, including the one-dimensional method, consistently perform among the best in the molecular test cases. Finally, the timing analysis shows the approximated one-dimensional method is far more efficient than standard second-order methods in the PBE force calculations. © 2017 Wiley Periodicals, Inc.

Modeling Membrane Protein-Ligand Binding Interactions: The Human Purinergic Platelet Receptor

Membrane proteins, due to their roles as cell receptors and signaling mediators, make prime candidates for drug targets. The computational analysis of protein-ligand binding affinities has been widely employed as a tool in rational drug design efforts. Although efficient implicit solvent-based methods for modeling globular protein-ligand binding have been around for many years, the extension of such methods to membrane protein-ligand binding is still in its infancy. In this study, we extended the widely used Amber/MMPBSA method to model membrane protein-ligand systems, and we used it to analyze protein-ligand binding for the human purinergic platelet receptor (P2Y12R), a prominent drug target in the inhibition of platelet aggregation for the prevention of myocardial infarction and stroke. The binding affinities, computed by the Amber/MMPBSA method using standard parameters, correlate well with experiment. A detailed investigation of these parameters was conducted to assess their impact on the accuracy of the method. These analyses show the importance of properly treating the nonpolar solvation interactions and the electrostatic polarization in the binding of nucleotide agonists and non-nucleotide antagonists to P2Y12R. On the basis of the crystal structures and the experimental conditions in the binding assay, we further hypothesized that the nucleotide agonists lose their bound magnesium ion upon binding to P2Y12R, and our computational study supports this hypothesis. Ultimately, this work illustrates the value of computational analysis in the interpretation of experimental binding reactions.

Calculating protein-ligand binding affinities with MMPBSA: Method and error analysis

Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) methods have become widely adopted in estimating protein-ligand binding affinities due to their efficiency and high correlation with experiment. Here different computational alternatives were investigated to assess their impact to the agreement of MMPBSA calculations with experiment. Seven receptor families with both high-quality crystal structures and binding affinities were selected. First the performance of nonpolar solvation models was studied and it was found that the modern approach that separately models hydrophobic and dispersion interactions dramatically reduces RMSD's of computed relative binding affinities. The numerical setup of the Poisson-Boltzmann methods was analyzed next. The data shows that the impact of grid spacing to the quality of MMPBSA calculations is small: the numerical error at the grid spacing of 0.5 Å is already small enough to be negligible. The impact of different atomic radius sets and different molecular surface definitions was further analyzed and weak influences were found on the agreement with experiment. The influence of solute dielectric constant was also analyzed: a higher dielectric constant generally improves the overall agreement with experiment, especially for highly charged binding pockets. The data also showed that the converged simulations caused slight reduction in the agreement with experiment. Finally the direction of estimating absolute binding free energies was briefly explored. Upon correction of the binding-induced rearrangement free energy and the binding entropy lost, the errors in absolute binding affinities were also reduced dramatically when the modern nonpolar solvent model was used, although further developments were apparently necessary to further improve the MMPBSA methods. © 2016 Wiley Periodicals, Inc.

Characterization of Promiscuous Binding of Phosphor Ligands to Breast-Cancer-Gene 1 (BRCA1) C-Terminal (BRCT): Molecular Dynamics, Free Energy, Entropy and Inhibitor Design

Inhibition of the protein-protein interaction (PPI) mediated by breast-cancer-gene 1 C-terminal (BRCT) is an attractive strategy to sensitize breast and ovarian cancers to chemotherapeutic agents that induce DNA damage. Such inhibitors could also be used for studies to understand the role of this PPI in DNA damage response. However, design of BRCT inhibitors is challenging because of the inherent flexibility associated with this domain. Several studies identified short phosphopeptides as tight BRCT binders. Here we investigated the thermodynamic properties of 18 phosphopeptides or peptide with phosphate mimic and three compounds with phosphate groups binding to BRCT to understand promiscuous molecular recognition and guide inhibitor design. We performed molecular dynamics (MD) simulations to investigate the interactions between inhibitors and BRCT and their dynamic behavior in the free and bound states. MD simulations revealed the key role of loops in altering the shape and size of the binding site to fit various ligands. The mining minima (M2) method was used for calculating binding free energy to explore the driving forces and the fine balance between configuration entropy loss and enthalpy gain. We designed a rigidified ligand, which showed unfavorable experimental binding affinity due to weakened enthalpy. This was because it lacked the ability to rearrange itself upon binding. Investigation of another phosphate group containing compound, C1, suggested that the entropy loss can be reduced by preventing significant narrowing of the energy well and introducing multiple new compound conformations in the bound states. From our computations, we designed an analog of C1 that introduced new intermolecular interactions to strengthen attractions while maintaining small entropic penalty. This study shows that flexible compounds do not always encounter larger entropy penalty, compared with other more rigid binders, and highlights a new strategy for inhibitor design.

Promiscuous proteins are commonly observed in biological systems, such as modular domains that recognize phosphopeptides during signal transduction. The use of phosphopeptides and compounds with phosphate groups as inhibitors to protein–protein interactions have attracted increasing interest for years. By using atomistic molecular dynamics simulations, we are able to perform detailed analyses of the dihedral space to explore protein fluctuation upon ligand binding to better understand promiscuous molecular recognition. Free energy calculation can further provide insights into the mechanism of binding, including both enthalpic and entropic contributions for molecular recognition, which assist in inhibitor design. Our calculation results show that pre-rigidifying a ligand is not always advantageous, suggesting the challenge in retaining optimized intermolecular interactions in pre-rigidified ligand. Instead, certain flexible ligands with multiple binding conformations can reduce entropic penalty, and therefore improves binding affinity. According to our computations, we can introduce new intermolecular interactions to flexible ligand to strengthen attractions while maintaining small entropic penalty by retaining its plasticity in the bound conformation. The study might cast light on a new general strategy for designing inhibitors targeting promiscuous modular domains and protein–protein interactions.

Base Flipping by MTERF1 Can Accommodate Multiple Conformations and Occurs in a Stepwise Fashion

Human mitochondrial transcription termination occurs within the leu-tRNA gene and is mediated by the DNA binding protein MTERF1. The crystal structure of MTERF1 bound to the canonical termination sequence reveals a rare base flipping event that involves the eversion of three nucleotides. These nucleotides are stabilized by stacking interactions with three MTERF1 residues, which are essential not only for base flipping but also for termination activity. To further understand the mechanism of base flipping, we examined each of the individual stacking interactions in structural, energetic and functional detail. Individual substitutions of Arg162, Tyr288 and Phe243 have revealed unequal contributions to overall termination activity. Furthermore, our work identifies an important role for Phe322 in the base flipping mechanism and we demonstrate how Phe322 and Phe243 are important for coupling base flipping between the heavy and light strand DNA chains. We propose a stepwise model for the base flipping process that recapitulates our observations. Finally, we show that MTERF1 has the ability to accommodate alternate active conformations. The adaptability of base flipping has implications for MTERF1 function and for the putative function of MTERF1 at alternative binding sites in human mitochondria.

Applications of MMPBSA to Membrane Proteins I: Efficient Numerical Solutions of Periodic Poisson-Boltzmann Equation

Continuum solvent models have been widely used in biomolecular modeling applications. Recently much attention has been given to inclusion of implicit membranes into existing continuum Poisson-Boltzmann solvent models to extend their applications to membrane systems. Inclusion of an implicit membrane complicates numerical solutions of the underlining Poisson-Boltzmann equation due to the dielectric inhomogeneity on the boundary surfaces of a computation grid. This can be alleviated by the use of the periodic boundary condition, a common practice in electrostatic computations in particle simulations. The conjugate gradient and successive over-relaxation methods are relatively straightforward to be adapted to periodic calculations, but their convergence rates are quite low, limiting their applications to free energy simulations that require a large number of conformations to be processed. To accelerate convergence, the Incomplete Cholesky preconditioning and the geometric multigrid methods have been extended to incorporate periodicity for biomolecular applications. Impressive convergence behaviors were found as in the previous applications of these numerical methods to tested biomolecules and MMPBSA calculations.

Switches of hydrogen bonds during ligand-protein association processes determine binding kinetics

Revealing the processes of ligand-protein associations deepens our understanding of molecular recognition and binding kinetics. Hydrogen bonds (H-bonds) play a crucial role in optimizing ligand-protein interactions and ligand specificity. In addition to the formation of stable H-bonds in the final bound state, the formation of transient H-bonds during binding processes contributes binding kinetics that define a ligand as a fast or slow binder, which also affects drug action. However, the effect of forming the transient H-bonds on the kinetic properties is little understood. Guided by results from coarse-grained Brownian dynamics simulations, we used classical molecular dynamics simulations in an implicit solvent model and accelerated molecular dynamics simulations in explicit waters to show that the position and distribution of the H-bond donor or acceptor of a drug result in switching intermolecular and intramolecular H-bond pairs during ligand recognition processes. We studied two major types of HIV-1 protease ligands: a fast binder, xk263, and a slow binder, ritonavir. The slow association rate in ritonavir can be attributed to increased flexibility of ritonavir, which yields multistep transitions and stepwise entering patterns and the formation and breaking of complex H-bond pairs during the binding process. This model suggests the importance of conversions of spatiotemporal H-bonds during the association of ligands and proteins, which helps in designing inhibitors with preferred binding kinetics.

Biological applications of classical electrostatics methods

Continuum electrostatics modeling of solvation based on the Poisson–Boltzmann (PB) equation has gained wide acceptance in biomolecular applications such as energetic analysis and structural visualization. Successful application of the PB solvent models requires careful calibration of the solvation parameters. Extensive testing and validation is also important to ensure accuracy in their applications. Limitation in the continuum modeling of solvation is also a known issue in certain biomolecular applications. Growing interest in membrane systems has further spurred developmental efforts to allow inclusion of membrane in the PB solvent models. Despite their past successes due to careful parameterization, algorithm development and parallel implementation, there is still much to be done to improve their transferability from the small molecular systems upon which they were developed and validated to complex macromolecular systems as advances in technology continue to push forward, providing ever greater computational resources to researchers to study more interesting biological systems of higher complexity.

Numerical Poisson-Boltzmann Model for Continuum Membrane Systems

Membrane protein systems are important computational research topics due to their roles in rational drug design. In this study, we developed a continuum membrane model utilizing a level set formulation under the numerical Poisson-Boltzmann framework within the AMBER molecular mechanics suite for applications such as protein-ligand binding affinity and docking pose predictions. Two numerical solvers were adapted for periodic systems to alleviate possible edge effects. Validation on systems ranging from organic molecules to membrane proteins up to 200 residues, demonstrated good numerical properties. This lays foundations for sophisticated models with variable dielectric treatments and second-order accurate modeling of solvation interactions.

Exploring a charge-central strategy in the solution of Poisson's equation for biomolecular applications

Continuum solvent treatments based on the Poisson-Boltzmann equation have been widely accepted for energetic analysis of biomolecular systems. In these approaches, the molecular solute is treated as a low dielectric region and the solvent is treated as a high dielectric continuum. The existence of a sharp dielectric jump at the solute-solvent interface poses a challenge to model the solvation energetics accurately with such a simple mathematical model. In this study, we explored and evaluated a strategy based on the "induced surface charge" to eliminate the dielectric jump within the finite-difference discretization scheme. In addition to the use of the induced surface charges in solving the equation, the second-order accurate immersed interface method is also incorporated to discretize the equation. The resultant linear system is solved with the GMRES algorithm to explicitly impose the flux conservation condition across the solvent-solute interface. The new strategy was evaluated on both analytical and realistic biomolecular systems. The numerical tests demonstrate the feasibility of utilizing induced surface charge in the finite-difference solution of the Poisson-Boltzmann equation. The analysis data further show that the strategy is consistent with theory and the classical finite-difference method on the tested systems. Limitations of the current implementations and further improvements are also analyzed and discussed to fully bring out its potential of achieving higher numerical accuracy.

Mechanistic insights into phosphopeptide--BRCT domain association: preorganization, flexibility, and phosphate recognition

Promiscuous proteins are commonly observed in biological systems, for example, in modular domains that recognize phosphopeptides during signal transduction. This promiscuous recognition is of fundamental interest in chemistry and biology but is challenging when designing phosphopeptides in silico for cell biology studies. To investigate promiscuous recognition and binding processes of phosphopeptides and the modular domain, we selected a domain essential in breast cancer-the breast-cancer-associated protein 1 (BRCA1) C-terminal (BRCT) repeats as our model system. We performed molecular dynamics simulations and detailed analyses of the dihedral space to study protein fluctuation and conformational changes with phosphopeptide binding. We also studied the association processes of phosphorylated and unphosphorylated peptides using Brownian dynamics with a coarse-grained model. We found that the BRCT domain is preorganized for phosphopeptide binding but has a moderate arrangement of side chains to form complexes with various types of phosphopeptides. Phosphopeptide binding restricts the system motion in general, while the nonpolar phosphopeptide becomes more flexible in the bound state. Our analysis found that the BRCT domain utilizes different mechanisms, usually termed lock and key, induced-fit, and population-shift/conformational-selection models, to recognize peptides with different features. Brownian dynamics simulations revealed that the charged phosphate group may not always accelerate peptide association processes, but it helps the phosphopeptide orient into binding pockets accurately and stabilizes the complex. This work provides insights into molecular recognition in the promiscuous protein system.

Reducing grid-dependence in finite-difference Poisson-Boltzmann calculations

Grid dependence in numerical reaction field energies and solvation forces is a well-known limitation in the finite-difference Poisson-Boltzmann methods. In this study we have investigated several numerical strategies to overcome the limitation. Specifically, we have included trimer arc dots during analytical molecular surface generation to improve the convergence of numerical reaction field energies and solvation forces. We have also utilized the level set function to trace the molecular surface implicitly to simplify the numerical mapping of the grid-independent solvent excluded surface. We have further explored to combine the weighted harmonic averaging of boundary dielectrics with a charge-based approach to improve the convergence and stability of numerical reaction field energies and solvation forces. Our test data show that the convergence and stability in both numerical energies and forces can be improved significantly when the combined strategy is applied to either the Poisson equation or the full Poisson-Boltzmann equation.

Mechanism of PhosphoThreonine/Serine Recognition and Specificity for Modular Domains from All-atom Molecular Dynamics

Background

Phosphopeptide-binding domains mediate many vital cellular processes such as signal transduction and protein recognition. We studied three well-known domains important for signal transduction: BRCT repeats, WW domain and forkhead-associated (FHA) domain. The first two recognize both phosphothreonine (pThr) and phosphoserine (pSer) residues, but FHA has high specificity for pThr residues. Here we used molecular dynamics (MD) simulations to reveal how FHA exclusively chooses pThr and how BRCT and WW recognize both pThr/pSer. The work also investigated the energies and thermodynamic information of intermolecular interactions.

Results

Simulations carried out included wide-type and mutated systems. Through analysis of MD simulations, we found that the conserved His residue defines dual loops feature of the FHA domain, which creates a small cavity reserved for only the methyl group of pThr. These well-organized loop interactions directly response to the pThr binding selectivity, while single loop (the 2nd phosphobinding site of FHA) or in combination with α-helix (BRCT repeats) or β-sheet (WW domain) fail to differentiate pThr/pSer.

Conclusions

Understanding the domain pre-organizations constructed by conserved residues and the driving force of domain-phosphopeptide recognition provides structural insight into pThr specific binding, which also helps in engineering proteins and designing peptide inhibitors.