Generalized Born Solvation Research Articles

High-Scalability Parallelization of a Molecular Modeling Application: Performance and Productivity Comparison Between OpenMP and MPI Implementations

Important components of molecular modeling applications are estimation and minimization of the internal energy of a molecule. For macromolecules such as proteins and amino acids, energy estimation is performed using empirical equations known as force fields. Over the past several decades, much effort has been directed towards improving the accuracy of these equations, and the resulting increased accuracy has come at the expense of greater computational complexity. For example, the interactions between a protein and surrounding water molecules have been modeled with improved accuracy using the generalized Born solvation model, which increases the computational complexity to O (n3). Fortunately, many force-field calculations are amenable to parallel execution. This paper describes the steps that were required to transform the Born calculation from a serial program into a parallel program suitable for parallel execution in both the OpenMP and MPI environments. Measurements of the parallel performance on a symmetric multiprocessor reveal that the Born calculation scales well for up to 144 processors. In some cases the OpenMP implementation scales better than the MPI implementation, but in other cases the MPI implementation scales better than the OpenMP implementation. However, in all cases the OpenMP implementation performs better than the MPI implementation, and requires less programming effort as well.

Folding simulations with novel conformational search method

A novel scheme for fast conformational search has been developed by combining the replica exchange method (REM) with the generalized effective potential concept. The new method, referred to Q-REM [S. Jang et al. Phys. Rev. Lett. 91, 058305 (2003)], is expected to provide a useful alternative to the conventional REM for effective conformational sampling of complex systems. The authors have performed folding simulations of the Trp-cage miniprotein using Q-REM. All atom level simulations with generalized Born solvent access-area solvation model show that successful folding can be observed with much smaller number of replicas in Q-REM compared to the conventional REM. It can be concluded that the new method has potential to significantly improve sampling efficiency, allowing simulations of more challenging systems.

Influence of temperature, friction, and random forces on folding of the B‐domain of staphylococcal protein A: All‐atom molecular dynamics in implicit solvent

The influences of temperature, friction, and random forces on the folding of protein A have been analyzed. A series of all-atom molecular dynamics folding simulations with the Amber ff99 potential and Generalized Born solvation, starting from the fully extended chain, were carried out for temperatures from 300 to 500 K, using (a) the Berendsen thermostat (with no explicit friction or random forces) and (b) Langevin dynamics (with friction and stochastic forces explicitly present in the system). The simulation temperature influences the relative time scale of the major events on the folding pathways of protein A. At lower temperatures, helix 2 folds significantly later than helices 1 and 3. However, with increasing temperature, the folding time of helix 2 approaches the folding times of helices 1 and 3. At lower temperatures, the complete formation of secondary and tertiary structure is significantly separated in time whereas, at higher temperatures, they occur simultaneously. These results suggest that some earlier experimental and theoretical observations of folding events, e.g., the order of helix formation, could depend on the temperature used in those studies. Therefore, the differences in temperature used could be one of the reasons for the discrepancies among published experimental and computational studies of the folding of protein A. Friction and random forces do not change the folding pathway that was observed in the simulations with the Berendsen thermostat, but their explicit presence in the system extends the folding time of protein A.

Second derivatives in generalized Born theory

Generalized Born solvation models offer a popular method of including electrostatic aspects of solvation free energies within an analytical model that depends only upon atomic coordinates, charges, and dielectric radii. Here, we describe how second derivatives with respect to Cartesian coordinates can be computed in an efficient manner that can be distributed over multiple processors. This approach makes possible a variety of new methods of analysis for these implicit solvation models. We illustrate three of these methods here: the use of Newton-Raphson optimization to obtain precise minima in solution; normal mode analysis to compute solvation effects on the mechanical properties of DNA; and the calculation of configurational entropies in the MM/GBSA model. An implementation of these ideas, using the Amber generalized Born model, is available in the nucleic acid builder (NAB) code, and we present examples for proteins with up to 45,000 atoms. The code has been implemented for parallel computers using both the OpenMP and MPI environments, and good parallel scaling is seen with as many as 144 OpenMP processing threads or MPI processing tasks.

New-Generation Amber United-Atom Force Field

We have developed a new-generation Amber united-atom force field for simulations involving highly demanding conformational sampling such as protein folding and protein-protein binding. In the new united-atom force field, all hydrogens on aliphatic carbons in all amino acids are united with carbons except those on Calpha. Our choice of explicit representation of all protein backbone atoms aims at minimizing perturbation to protein backbone conformational distributions and to simplify development of backbone torsion terms. Tests with dipeptides and solvated proteins show that our goal is achieved quite successfully. The new united-atom force field uses the same new RESP charging scheme based on B3LYP/cc-pVTZ//HF/6-31g** quantum mechanical calculations in the PCM continuum solvent as that in the Duan et al. force field. van der Waals parameters are empirically refitted starting from published values with respect to experimental solvation free energies of amino acid side-chain analogues. The suitability of mixing new point charges and van der Waals parameters with existing Amber covalent terms is tested on alanine dipeptide and is found to be reasonable. Parameters for all new torsion terms are refitted based on the new point charges and the van der Waals parameters. Molecular dynamics simulations of three small globular proteins in the explicit TIP3P solvent are performed to test the overall stability and accuracy of the new united-atom force field. Good agreements between the united-atom force field and the Duan et al. all-atom force field for both backbone and side-chain conformations are observed. In addition, the per-step efficiency of the new united-atom force field is demonstrated for simulations in the implicit generalized Born solvent. A speedup around two is observed over the Duan et al. all-atom force field for the three tested small proteins. Finally, the efficiency gain of the new united-atom force field in conformational sampling is further demonstrated with a well-known toy protein folding system, an 18 residue polyalanine in distance-dependent dielectric. The new united-atom force field is at least a factor of 200 more efficient than the Duan et al. all-atom force field for ab initio folding of the tested peptide.

Is Poisson-Boltzmann theory insufficient for protein folding simulations?

The Poisson-Boltzmann theory has been widely used in the studies of energetics and conformations of biological macromolecules. Recently, introduction of the efficient generalized Born approximation has greatly extended its applicability to areas such as protein folding simulations where highly efficient computation is crucial. However, limitations have been found in the folding simulations of a well-studied beta hairpin with several generalized Born implementations and different force fields. These studies have raised the question whether the underlining Poisson-Boltzmann theory, on which the generalized Born model is calibrated, is adequate in the treatment of polar interactions for the challenging protein folding simulations. To address the question whether the Poisson-Boltzmann theory in the current formalism might be insufficient, we directly tested our efficient numerical Poisson-Boltzmann implementation in the beta-hairpin folding simulation. Good agreement between simulation and experiment was found for the beta-hairpin equilibrium structures when the numerical Poisson-Boltzmann solvent and a recently improved generalized Born solvent were used. In addition simulated thermodynamic properties also agree well with experiment in both solvents. Finally, an overall agreement on the beta-hairpin folding mechanism was found between the current and previous studies. Thus, our simulations indicate that previously observed limitations are most likely due to imperfect calibration in previous generalized Born models but not due to the limitation of the Poisson-Boltzmann theory.

Investigation of Salt Bridge Stability in a Generalized Born Solvent Model.

Potentials of mean force (PMFs) of salt bridge formation between oppositely charged amino acid side chains were calculated both in explicit solvent and in a Generalized Born (GB) continuum solvent model to quantify the potential overstabilization of side chain ion pairs in GB relative to explicit solvation. These show that salt bridges are too stable by as much as 3-4 kcal/mol in the GB solvent models that we tested, consistent with previously reported observations of significantly different structural ensembles in GB models and explicit solvent for proteins containing ionizable groups. We thus investigated a simple empirical correction, wherein the intrinsic GB radii of hydrogen atoms bound to charged nitrogen atoms are reduced, effectively increasing the desolvation penalty of the positively charged groups. The thermodynamics of salt bridge formation were considerably improved, as exemplified by the close match of the corrected GB PMF to the reference explicit solvent PMF, and more significantly by our ability to closely reproduce the experimental temperature melting profile of the TC5b Trp-cage miniprotein, which is otherwise highly distorted by prevalent non-native salt bridges when using standard GB parameters.

PKa of Acetate in Water: A Computational Study

Several computational methods including the conductor-like polarizable continuum model, CPCM with both UAKS and UAHF cavities, Cramer and Truhlar's generalized Born solvation model, SM5.4(AM1), SM5.4(PM3), and SM5.43R(mPW1PW91/6-31+G(d)), and mixed QM/MM-Ewald simulations were used to calculate the pK(a) values of acetate and bicarbonate anions in aqueous solution. This work provided a critical and comprehensive assessment of the quality of these theoretical models in the calculation of aqueous solvation free energies for the singly charged acetate and bicarbonate ions, as well as the doubly charged acetate dianion and carbonate dianion. It was shown that QM/MM-Ewald simulations could give an accurate and consistent evaluation of the pK(a) values of acetate and bicarbonate based on both the relative and absolute pK(a) formulas, while other methods could yield satisfactory results only for certain calculations. However, this does not mean that the current QM/MM-Ewald protocol is superior to other methods. The useful information obtained in this investigation is that both the absolute and relative pK(a) formulas should better be tested in accurate calculations of pK(a) values based on any methods.

Virtual Ligand Screening against Escherichia coli Dihydrofolate Reductase: Improving Docking Enrichment Using Physics-Based Methods

Motivated by their participation in the McMaster Data-Mining and Docking Competition, the authors developed 2 new computational technologies and applied them to docking against Escherichia coli dihydrofolate reductase: a receptor preparation procedure that incorporates rotamer optimization of side chains and a physics-based rescoring procedure for estimating relative binding affinities of the protein-ligand complexes. Both methods use the same energy function, consisting of the all-atom OPLS-AA force field and a generalized Born solvent model, which treats the protein receptor and small-molecule ligands in a consistent manner. Thus, the energy function is similar to that used in more sophisticated approaches, such as free-energy perturbation and the molecular mechanics Poisson-Boltzmann/surface area, but sampling during the rescoring procedure is limited to simple energy minimization of the ligand. The use of a highly efficient minimization algorithm permitted the authors to apply this rescoring procedure to hundreds of thousands of protein-ligand complexes during the competition, using a modest Linux cluster. To test these methods, they used the 12 competitive inhibitors identified in the training set, plus methotrexate, as positive controls in enrichment studies with both the training and test sets, each containing 50,000 compounds. The key conclusion is that combining the receptor preparation and rescoring methods makes it possible to identify most of the positive controls within the top few tenths of a percent of the rank-ordered training and test set libraries.

Exploring Phase‐Transfer Catalysis with Molecular Dynamics and 3D/4D Quantitative Structure—Selectivity Relationships.

Quantitative Structure-Selectivity Relationships (QSSR) are developed for a library of 40 phase-transfer asymmetric catalysts, based around quaternary ammonium salts, using Comparative Molecular Field Analysis (CoMFA) and closely related variants. Due to the flexibility of these catalysts, we use molecular dynamics (MD) with an implicit Generalized Born solvent model to explore their conformational space. Comparison with crystal data indicates that relevant conformations are obtained and that, furthermore, the correct biphenyl twist conformation is predicted, as illustrated by the superiority of the resulting model (leave-one-out q(2) = 0.78) compared to a random choice of low-energy conformations for each catalyst (average q(2) = 0.22). We extend this model by incorporating the MD trajectory directly into a 4D QSSR and by Boltzmann-weighting the contribution of selected minimized conformations, which we refer to as '3.5D' QSSR. The latter method improves on the predictive ability of the 3D QSSR (leave-one-out q(2) = 0.83), as confirmed by repeated training/test splits.

Exploring Phase-Transfer Catalysis with Molecular Dynamics and 3D/4D Quantitative Structure−Selectivity Relationships

Quantitative Structure-Selectivity Relationships (QSSR) are developed for a library of 40 phase-transfer asymmetric catalysts, based around quaternary ammonium salts, using Comparative Molecular Field Analysis (CoMFA) and closely related variants. Due to the flexibility of these catalysts, we use molecular dynamics (MD) with an implicit Generalized Born solvent model to explore their conformational space. Comparison with crystal data indicates that relevant conformations are obtained and that, furthermore, the correct biphenyl twist conformation is predicted, as illustrated by the superiority of the resulting model (leave-one-out q(2) = 0.78) compared to a random choice of low-energy conformations for each catalyst (average q(2) = 0.22). We extend this model by incorporating the MD trajectory directly into a 4D QSSR and by Boltzmann-weighting the contribution of selected minimized conformations, which we refer to as '3.5D' QSSR. The latter method improves on the predictive ability of the 3D QSSR (leave-one-out q(2) = 0.83), as confirmed by repeated training/test splits.

Misfolded free energy surface of a peptide with αββ motif (1PSV) using the generalized Born solvation model

Using the replica exchange molecular dynamics method, we carried out ab initio folding study of 1PSV using AMBER param99MOD2 and the generalized Born (GB) solvation model. A total of 10 ns run for each replica resulted in a distorted free energy surface populated by two major misfolded conformers. The native structure of 1PSV was never located in the free energy surface. This distortion is mainly due to the defect of the current energy model including the GB solvation term. Therefore, the current energy model should be further improved for correct folding predictions of the proteins with a αββ motif.

Simulation of the Interaction Between ScyTx and Small Conductance Calcium-Activated Potassium Channel by Docking and MM-PBSA

Computational methods are employed to simulate interaction of scorpion toxin ScyTx in complex with the small conductance calcium-activated potassium channel rsk2. All of available 25 structures of ScyTx in the Protein Data Bank determined by NMR were considered for improving performance of rigid protein docking of ZDOCK. Four main binding modes were found among a large number of predicted complexes by using clustering analysis, screening with expert knowledge, energy minimization, and molecular dynamics simulations. The quality and validity of the resulting complexes were further evaluated by molecular dynamics simulations with the generalized Born solvation model and by calculation of relative binding free energies with the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) in the AMBER 7 suit of programs. The complex formed by the 22nd structure of the ScyTx and rsk2 channel was identified as the most favorable complex by using a combination of computational methods, which contain further introduction of flexibility without restraining residue side chain. From the resulted spatial structure of the ScyTx and rsk2 channel, ScyTx associates the mouth of the rsk2 channel with α-helix rather than β-sheet. Structural analysis first revealed that Arg 13 played a novel and vital role of blocking the pore of the rsk2 channel, whose role is remarkably different from that of highly homologous scorpion toxin P05. Between the interfaces in the ScyTx-rsk2 complex, strong electrostatic interaction and hydrogen bonds exist between Arg 13 of ScyTx and Gly-Tyr-Gly-Asp sequential residues located in the four symmetrical chains of the pore region. Simultaneously, five hydrogen bonds between Arg 6 of ScyTx and Asp 341(C), Val 366(C), and Pro 367(C), and electrostatic interaction between Arg 6 of ScyTx and Asp 364(B) and Asp 341(C) are also found by structural analysis. In addition, His 31 located at the C-terminal of ScyTx is surrounded by Val 342(A), Asp 364(A), Met 365(A), Pro 367(B), and Asn 366(B) within a contact distance of 4.0 Å. These simulation results are in good agreement with experimental data and can effectively explain the binding phenomena between ScyTx and the potassium channel at the level of molecular spatial structure. The consistency between results of molecular modeling and experimental data strongly suggests that our spatial structure model of the ScyTx-rsk2 complex is reasonable. Therefore, molecular docking combined with molecular dynamics simulations followed by molecular mechanics Poisson-Boltzmann surface area analysis is an attractive approach for modeling scorpion toxin-potassium channel complexes a priori for further biological studies.

Molecular Dynamics Simulations of DNA Using the Generalized Born Solvation Model: Quantitative Comparisons with Explicit Solvation Results

It is now well established that, with the correct treatment of hydration and long-range electrostatic effects, simulations of nucleic acid structure, dynamics, and recognition can yield remarkably reliable results and have genuine predictive power. While the “gold standard” remains modeling studies with the explicit treatment of water molecules and counterions, new implicit hydration models are being advanced as suitable choices for simulations where a shorter time-to-solution is required. In this study, we have evaluated the reliability of one such model, the generalized Born/surface area (GB/SA) continuum model by using it to study the interaction between the minor groove binding drug Hoechst 33258 and the DNA dodecamer d(CTTTTGCAAAAG)2. There is extensive NMR and explicit hydrated modeling data on this system, and it shows the interesting property of highly cooperative binding of two molecules of the drug, such that the 1:1 complex is never observed. We find that with a suitable choice of parameters, the GB/SA model can perform adequately, reproducing the structural and dynamical characteristics observed in the previous explicitly solvated simulations in approximately a quarter of the computational time. However, limitations of this GB/SA method become apparent when the thermodynamic properties are evaluated.

A Free Energy Based Computational Pathway from Chemical Templates to Lead Compounds: A Case Study of COX-2 Inhibitors

Automation of lead compound design in silico given the structure of the protein target and a definition of its active site vies for the top of the wish list in any drug discovery programme. We present here an enumeration of steps starting from chemical templates and propose a solution at the state of the art, in the form of a system independent comprehensive computational pathway. This methodology is illustrated with cyclooxygenase-2 (COX-2) as a target. We built candidate molecules including a few Non Steroidal Anti-inflammatory Drugs (NSAIDs) from chemical templates, passed them through empirical filters to assess drug-like properties, optimized their geometries, derived partial atomic charges via quantum calculations, performed Monte Carlo docking, carried out molecular mechanics and developed free energy estimates with Molecular Mechanics Generalized Born Solvent Accessibility (MMGBSA) methodology for each of the candidate molecules. For the case of aspirin, we also conducted molecular dynamics on the enzyme, the drug and the complex with explicit solvent followed by binding free energy analysis. Collectively, the results obtained from the above studies viz. sorting of drugs from non-drugs, semi-quantitative estimates of binding free energies, amply demonstrate the viability of the strategy proposed for lead selection/design for biomolecular targets.

Molecular dynamics and free energy analysis of neuraminidase-ligand interactions.

We report molecular dynamics calculations of neuraminidase in complex with an inhibitor, 4-amino-2-deoxy-2,3-didehydro-N-acetylneuraminic acid (N-DANA), with subsequent free energy analysis of binding by using a combined molecular mechanics/continuum solvent model approach. A dynamical model of the complex containing an ionized Glu119 amino acid residue is found to be consistent with experimental data. Computational analysis indicates a major van der Waals component to the inhibitor-neuraminidase binding free energy. Based on the N-DANA/neuraminidase molecular dynamics trajectory, a perturbation methodology was used to predict the binding affinity of related neuraminidase inhibitors by using a force field/Poisson-Boltzmann potential. This approach, incorporating conformational search/local minimization schemes with distance-dependent dielectric or generalized Born solvent models, correctly identifies the most potent neuraminidase inhibitor. Mutation of the key ligand four-substituent to a hydrogen atom indicates no favorable binding free energy contribution of a hydroxyl group; conversely, cationic substituents form favorable electrostatic interactions with neuraminidase. Prospects for further development of the method as an analysis and rational design tool are discussed.

NMR structure of the (+)-CPI-indole/d(GACTAATTGAC)-d(GTCAATTAGTC) covalent complex.

We report the NMR solution structure of (+)-CPI-indole (CPI, 1,2,8,8a-tetrahydrocyclopropa[c]pyrrolo[3,2-e]indol-4(5H)-one), an agent belonging to the CC-1065/duocarmycin family of antitumor compounds. This (+)-CPI-indole structure is covalently bound to d(G(1)ACTAATTGTC(11))-d(G(12)TCAATTAGTC(22)), a synthetic DNA duplex containing a high-affinity binding site. The three-dimensional structure has been determined by several cycles of restrained molecular dynamics calculations with a total of 563 NMR-derived constraints, both in vacuo and by using the generalized Born solvent continuum model. In-depth analysis of the structure of this ligand-DNA complex led to a detailed knowledge of the bound state conformation of the CPI-indole, the most simplified agent related to CC-1065 and duocarmycins, the parent members of a family of extremely potent antitumor compounds. Comparison of the CPI-indole bound conformation with those previously found for (+)-duocarmycin SA (DSA), its unnatural enantiomer (-)-DSA, and the demethoxylated analogue (+)-DSI in their DNA complexes provided additional evidence of the tight correlation between the catalytic effect exerted by DNA on the alkylation reaction and the extent of angular twist between the two planar heteroaromatic subunits of these agents. Additionally, comparison of the structural features of the DNA-bound state of a "naked" ligand, such as CPI-indole, with those of various other duocarmycin agents provided useful information for the interpretation of the observed effects on chemical reactivity of the different substitution patterns at the hemispheres of these types of complex.

Stability of the I-motif Structure Is Related to the Interactions between Phosphodiester Backbones

The i-motif DNA tetrameric structure is formed of two parallel duplexes intercalated in a head-to-tail orientation, and held together by hemiprotonated cytosine pairs. The four phosphodiester backbones forming the structure define two narrow and wide grooves. The short interphosphate distances across the narrow groove induce a strong repulsion which should destabilize the tetramer. To investigate this point, molecular dynamics simulations were run on the [d(C2)]4 and [d(C4)]4 tetramers in 3′E and 5′E topologies, for which the interaction of the phosphodiester backbones through the narrow groove is different. The analysis of the simulations, using the Molecular Mechanics Generalized Born Solvation Area and Molecular Mechanics Poisson-Boltzmann Solvation Area approaches, shows that it is the van der Waals energy contribution which displays the largest relative difference between the two topologies. The comparison of the solvent-accessible area of each topology reveals that the sugar-sugar interactions account for the greater stability of the 3′E topology. This stresses the importance of the sugar-sugar contacts across the narrow groove which, enforcing the optimal backbone twisting, are essential to the base stacking and the i-motif stability. Tighter interactions between the sugars are observed in the case of N-type sugar puckers.

Solution Structure of the Chick TGFβ Type II Receptor Ligand-binding Domain

The transforming growth factor β (TGFβ) signaling pathway influences cell proliferation, immune responses, and extracellular matrix reorganization throughout the vertebrate life cycle. The signaling cascade is initiated by ligand-binding to its cognate type II receptor. Here, we present the structure of the chick type II TGFβ receptor determined by solution NMR methods. Distance and angular constraints were derived from 15N and 13C edited NMR experiments. Torsion angle dynamics was used throughout the structure calculations and refinement. The 20 final structures were energy minimized using the generalized Born solvent model. For these 20 structures, the average backbone root-mean-square distance from the average structure is below 0.6 Å. The overall fold of this 109-residue domain is conserved within the superfamily of these receptors. Chick receptors fully recognize and respond to human TGFβ ligands despite only 60% identity at the sequence level. Comparison with the human TGFβ receptor determined by X-ray crystallography reveals different conformations in several regions. Sequence divergence and crystal packing interactions under low pH conditions are likely causes. This solution structure identifies regions were structural changes, however subtle, may occur upon ligand-binding. We also identified two very well conserved molecular surfaces. One was found to bind ligand in the crystallized human TGFβ3:TGFβ type II receptor complex. The other, newly identified area can be the interaction site with type I and/or type III receptors of the TGFβ signaling complex.

Prediction of helical peptide folding in an implicit water by a new molecular dynamics scheme with generalized effective potential

Previously we demonstrated the effectiveness of the recently developed q-jumping molecular dynamics simulation method (q-jumping MD) in vacuo for fast conformational searches and optimization purposes. In this work we attempt to further investigate the conformational searching capability of this new scheme by applying it to folding problems of helical peptides containing fully charged side chains in an implicit water. With a slightly modified q-jumping MD method using the all-atom empirical force field and its generalized Born solvation model, the current simulations at T=300 K all lead to fast helix folding with broad potential energy fluctuations, starting from their extended (linear) conformations. The present study demonstrates that this new MD scheme greatly enhances the rate of conformational changes, making it possible to explore low energy conformations of peptides in aqueous environments in a reasonably short time. Therefore, the all-atom based theoretical prediction of native solution structures of more challenging systems, such as helix bundles, β-sheets, and even small proteins may be a realistic possibility.