Enzyme-inhibitor association thermodynamics: explicit and continuum solvent studies

Adaptive Multilevel Splitting Method for Molecular Dynamics Calculation of Benzamidine-Trypsin Dissociation Time

Adaptive multilevel splitting (AMS) is a rare event sampling method that requires minimal parameter tuning and allows unbiased sampling of transition pathways of a given rare event. Previous simulation studies have verified the efficiency and accuracy of AMS in the calculation of transition times for simple systems in both Monte Carlo and molecular dynamics (MD) simulations. Now, AMS is applied for the first time to an MD simulation of protein-ligand dissociation, representing a leap in complexity from the previous test cases. Of interest is the dissociation rate, which is typically too low to be accessible to conventional MD. The present study joins other recent efforts to develop advanced sampling techniques in MD to calculate dissociation rates, which are gaining importance in the pharmaceutical field as indicators of drug efficacy. The system investigated here, benzamidine bound to trypsin, is an example common to many of these efforts. The AMS estimate of the dissociation rate was found to be (2.6 ± 2.4) × 10(2) s(-1), which compares well with the experimental value.

Rapid estimation of hydration thermodynamics of macromolecular regions

This work describes a novel protocol to efficiently calculate the local free energy of hydration of specific regions in macromolecules. The method employs Monte Carlo simulations in the grand canonical ensemble to generate water configurations in a selected spherical region in the macromolecule. Excess energy and entropy of hydration are calculated by analyzing the water configurational distributions following the recently published grid inhomogeneous solvation theory method [C. N. Nguyen, T. K. Young, and M. K. Gilson, J. Chem. Phys. 137, 044101 (2012)]. Our method involves the approximations of treating the macromolecule and distant solvent as rigid and performing calculations on multiple such conformations to account for conformational diversity. These approximations are tested against water configurations obtained from a molecular dynamics simulation. The method is validated by predicting the number and location of water molecules in 5 pockets in the protein Interleukin-1β for which experimental water occupancy data are available. Free energy values are validated against decoupling free energy perturbation calculations. The results indicate that the approximations used in the method enable efficient prediction of free energies of water displacement.

Placevent: an algorithm for prediction of explicit solvent atom distribution-application to HIV-1 protease and F-ATP synthase

We have created a simple algorithm for automatically predicting the explicit solvent atom distribution of biomolecules. The explicit distribution is coerced from the three-dimensional (3D) continuous distribution resulting from a 3D reference interaction site model (3D-RISM) calculation. This procedure predicts optimal location of solvent molecules and ions given a rigid biomolecular structure and the solvent composition. We show examples of predicting water molecules near the KNI-272 bound form of HIV-1 protease and predicting both sodium ions and water molecules near the rotor ring of F-adenosine triphosphate (ATP) synthase. Our results give excellent agreement with experimental structure with an average prediction error of 0.39-0.65 Å. Further, unlike experimental methods, this method does not suffer from the partial occupancy limit. Our method can be performed directly on 3D-RISM output within minutes. It is extremely useful for examining multiple specific solvent-solute interactions, as a convenient method for generating initial solvent structures for molecular dynamics calculations, and may assist in refinement of experimental structures. © 2012 Wiley Periodicals, Inc.

Complete reconstruction of an enzyme-inhibitor binding process by molecular dynamics simulations

The understanding of protein-ligand binding is of critical importance for biomedical research, yet the process itself has been very difficult to study because of its intrinsically dynamic character. Here, we have been able to quantitatively reconstruct the complete binding process of the enzyme-inhibitor complex trypsin-benzamidine by performing 495 molecular dynamics simulations of free ligand binding of 100 ns each, 187 of which produced binding events with an rmsd less than 2 Å compared to the crystal structure. The binding paths obtained are able to capture the kinetic pathway of the inhibitor diffusing from solvent (S0) to the bound (S4) state passing through two metastable intermediate states S2 and S3. Rather than directly entering the binding pocket the inhibitor appears to roll on the surface of the protein in its transition between S3 and the final binding pocket, whereas the transition between S2 and the bound pose requires rediffusion to S3. An estimation of the standard free energy of binding gives ΔG° = -5.2 ± 0.4 kcal/mol (cf. the experimental value -6.2 kcal/mol), and a two-states kinetic model k(on) = (1.5 ± 0.2) × 10(8) M(-1) s(-1) and k(off) = (9.5 ± 3.3) × 10(4) s(-1) for unbound to bound transitions. The ability to reconstruct by simple diffusion the binding pathway of an enzyme-inhibitor binding process demonstrates the predictive power of unconventional high-throughput molecular simulations. Moreover, the methodology is directly applicable to other molecular systems and thus of general interest in biomedical and pharmaceutical research.

Trypsin-ligand binding free energies from explicit and implicit solvent simulations with polarizable potential

We have calculated the binding free energies of a series of benzamidine-like inhibitors to trypsin with a polarizable force field using both explicit and implicit solvent approaches. Free energy perturbation has been performed for the ligands in bulk water and in protein complex with molecular dynamics simulations. The binding free energies calculated from explicit solvent simulations are well within the accuracy of experimental measurement and the direction of change is predicted correctly in all cases. We analyzed the molecular dipole moments of the ligands in gas, water and protein environments. Neither binding affinity nor ligand solvation free energy in bulk water shows much dependence on the molecular dipole moments of the ligands. Substitution of the aromatic or the charged group in the ligand results in considerable change in the solvation energy in bulk water and protein whereas the binding affinity varies insignificantly due to cancellation. The effect of chemical modification on ligand charge distribution is mostly local. Replacing benzene with diazine has minimal impact on the atomic multipoles at the amidinium group. We have also utilized an implicit solvent based end-state approach to evaluate the binding free energies of these inhibitors. In this approach, the polarizable multipole model combined with Poisson-Boltzmann/surface area (PMPB/SA) provides the electrostatic interaction energy and the polar solvation free energy. Overall the relative binding free energies obtained from the MM-PMPB/SA model are in good agreement with the experimental data.

Computational simulations of the Trichoderma reesei cellobiohydrolase I acting on microcrystalline cellulose Ibeta: the enzyme-substrate complex

Cellobiohydrolases are the dominant components of the commercially relevant Trichoderma reesei cellulase system. Although natural cellulases can totally hydrolyze crystalline cellulose to soluble sugars, the current enzyme loadings and long digestion times required render these enzymes less than cost effective for biomass conversion processes. It is clear that cellobiohydrolases must be improved via protein engineering to reduce processing costs. To better understand cellobiohydrolase function, new simulations have been conducted using charmm of cellobiohydrolase I (CBH I) from T.reesei interacting with a model segment (cellodextrin) of a cellulose microfibril in which one chain from the substrate has been placed into the active site tunnel mimicking the hypothesized configuration prior to final substrate docking (i.e., the +1 and +2 sites are unoccupied), which is also the structure following a catalytic bond scission. No tendency was found for the protein to dissociate from or translate along the substrate surface during this initial simulation, nor to align with the direction of the cellulose chains. However, a tendency for the decrystallized cellodextrin to partially re-anneal into the cellulose surface hints that the arbitrary starting configuration selected was not ideal.

Standard Free Energy of Binding from a One-Dimensional Potential of Mean Force

A practical approach that enables one to calculate the standard free energy of binding from a one-dimensional potential of mean force (PMF) is proposed. Umbrella sampling and the weighted histogram analysis method are used to generate a PMF along the reaction coordinate of binding. At each point, a restraint is applied orthogonal to the reaction coordinate to make possible the determination of the volume sampled by the ligand. The free energy of binding from an arbitrary unbound volume to the restrained bound form is calculated from the ratio of the PMF integrated over the bound region to that of the unbound. Adding the free energy changes from the standard-state volume to the unbound volume and from the restrained to the unrestrained bound state gives the standard free energy of binding. Exploration of the best choice of binding paths is also made. This approach is first demonstrated on a model binding system and then tested on the benzamidine-trypsin system for which reasonable agreement with experiment is found. A comparison is made with other methods to obtain the standard free energy of binding from the PMF.

Calculation of absolute ligand binding free energy to a ribosome-targeting protein as a function of solvent model

A comparative analysis is provided of the effect of different solvent models on the calculation of a potential of mean force (PMF) for determining the absolute binding affinity of the small molecule inhibitor pteroic acid bound to ricin toxin A-chain (RTA). Solvent models include the distance-dependent dielectric constant, several different generalized Born (GB) approximations, and a hybrid explicit/GB-based implicit solvent model. We found that the simpler approximation of dielectric screening and a GB model, with Born radii fitted to a switching-window dielectric-boundary surface Poisson solvent model, severely overpredicted the binding affinity as compared to the experimental value, estimated to range from -4.4 to -6.0 kcal/mol. In contrast, GB models that are parametrized to fit the Lee-Richards molecular surface performed much better, predicting binding free energy within 1-3 kcal/mol of experimental estimates. However, the predicted free-energy profiles of these GB models displayed alternative binding modes not observed in the crystal structure. Finally, the most rigorous and computationally costly approach in this work, which used a hybrid explicit/implicit solvent model, correctly determined a binding funnel in the PMF near the crystallographic bound state and predicted an absolute binding affinity that was 2 kcal/mol more favorable than the estimated experimental binding affinity.

"Tetrahedrality" and the relationship between collective structure and radial distribution functions in liquid water

Molecular dynamics simulations of pure liquid water under ambient conditions using four common empirical water models have been analyzed to determine how well the oxygen-oxygen radial distribution function, g(OO)(r), used as the sole criterion of congruence with experiment, captures variations in the actual anisotropic collective structuring for these models. The largest systematic deviations from tetrahedrality were found to be due to deformations of the angle between the two closet hydrogen bond donor neighbors, but for intrinsic geometric reasons, these were found to contribute less to g(OO)(r) than deformations of the angles between one hydrogen bond donor neighbor and one hydrogen bond acceptor neighbor. Relying exclusively on a qualitative characterization of the second peak in g(OO)(r) seems to overemphasize the differences between the structuring in some of these models.

Odorant binding and conformational dynamics in the odorant-binding protein

In mammals, the olfactory epithelium secretes odorant-binding proteins (OBPs), which are lipocalins found freely dissolved in the mucus layer protecting the olfactory neurons. OBPs may act as passive transporters of predominantly hydrophobic odorant molecules across the aqueous mucus layer, or they may play a more active role in which the olfactory neuronal receptor recognizes the OBP-ligand complex. To better understand the molecular events accompanying the initial steps in the olfaction process, we have performed molecular dynamics studies of rat and pig OBPs with the odorant molecule thymol. These calculations provide an atomic level description of conformational changes and pathway intermediates that remain difficult to study directly. A series of eight independent molecular dynamics trajectories of rat OBP permitted the observation of a consensus pathway for ligand unbinding and the calculation of the potential of mean force (PMF) along this path. Titration microcalorimetry confirmed the specific binding of thymol to this protein with a strong hydrophobic component. In both rat and pig OBPs we observed lipocalin strand pair opening in the presence of ligand, consistent with potential roles of these proteins in olfactive receptor recognition.

Calculation of absolute protein-ligand binding affinity using path and endpoint approaches

A comparative analysis is provided of rigorous and approximate methods for calculating absolute binding affinities of two protein-ligand complexes: the FKBP protein bound with small molecules 4-hydroxy-2-butanone and FK506. Our rigorous approach is an umbrella sampling technique where a potential of mean force is determined by pulling the ligand out of the protein active site over several simulation windows. The results of this approach agree well with experimentally observed binding affinities. Also assessed is a commonly used approximate endpoint approach, which separately estimates enthalpy, solvation free energy, and entropy. We show that this endpoint approach has numerous variations, all of which are prone to critical shortcomings. For example, conventional harmonic and quasiharmonic entropy estimation procedures produce disparate results for the relatively simple protein-ligand systems studied in this work.

Computer simulation study of intermolecular voids in unsaturated phosphatidylcholine lipid bilayers

Computer simulation of the liquid crystalline phase of five different hydrated unsaturated phosphadidylcholine (PC) lipid bilayers, i.e., membranes built up by 18:0/18:1omega9cis PC, 18:0/18:2omega6cis PC, 18:0/18:3omega3cis PC, 18:0/20:4omega6cis PC, and 18:0/22:6omega3cis PC molecules have been performed on the isothermal-isobaric ensemble at 1 atm and 303 K. (The notation n:domegapcis specifies the lipid tails: n refers to the total number of carbon atoms in the chain, d is the number of the methylene-interrupted double bonds, p denotes the number of carbons between the chain terminal CH(3) group and the nearest double bond, and cis refers to the conformation around the double bonds.) The characteristics of the free volume in these systems have been analyzed by means of a generalized version of the Voronoi-Delaunay method [M. G. Alinchenko et al., J. Phys. Chem. B 108, 19056 (2004)]. As a reference system, the hydrated bilayer of the saturated 14:014:0 PC molecules (dimyristoylphosphatidylcholine) has also been analyzed. It has been found that the profiles of the fraction of the free volume across the membrane exhibit a rather complex pattern. This fine structure of the free volume fraction profiles can be interpreted by dividing the membrane into three separate major zones (i.e., zones of the aqueous, polar, and apolar parts of the membrane) and defining five subzones within these zones according to the average position of various atomic groups in the membrane. The fraction of the free volume in the middle of the membrane is found to increase with increasing unsaturation of the sn-2 chain of the lipid molecule. This is due to the fact that with increasing number of methylene-interrupted double bonds the lipid tails become more flexible, and hence they do not extend to the middle of the membrane. It is found that there are no broad enough preformed channels in the bilayers through which small penetrants, such as water molecules, can readily go through; however, the existing channels can largely facilitate the permeation of these molecules.

Theoretical identification of proton channels in the quinol oxidase aa3 from Acidianus ambivalens

Heme-copper oxidases are membrane proteins found in the respiratory chain of aerobic organisms. They are the terminal electron acceptors coupling the translocation of protons across the membrane with the reduction of oxygen to water. Because the catalytic process occurs in the heme cofactors positioned well inside the protein matrix, proton channels must exist. However, due to the high structural divergence among this kind of proteins, the proton channels previously described are not necessarily conserved. In this work we modeled the structure of the quinol oxidase from Acidianus ambivalens using comparative modeling techniques for identifying proton channels. Additionally, given the high importance that water molecules may have in this process, we have developed a methodology, within the context of comparative modeling, to identify high water probability zones and to deconvolute them into chains of ordered water molecules. From our results, and from the existent information from other proteins from the same superfamily, we were able to suggest three possible proton channels: one K-, one D-, and one Q-spatial homologous proton channels. This methodology can be applied to other systems where water molecules are important for their biological function.

Molecular dynamics simulations of peptides and proteins with a continuum electrostatic model based on screened Coulomb potentials

A continuum electrostatics approach for molecular dynamics (MD) simulations of macromolecules is presented and analyzed for its performance on a peptide and a globular protein. The approach incorporates the screened Coulomb potential (SCP) continuum model of electrostatics, which was reported earlier. The model was validated in a broad set of tests some of which were based on Monte Carlo simulations that included single amino acids, peptides, and proteins. The implementation for large-scale MD simulations presented in this article is based on a pairwise potential that makes the electrostatic model suitable for fast analytical calculation of forces. To assess the suitability of the approach, a preliminary validation is conducted, which consists of (i) a 3-ns MD simulation of the immunoglobulin-binding domain of streptococcal protein G, a 56-residue globular protein and (ii) a 3-ns simulation of Dynorphin, a biological peptide of 17 amino acids. In both cases, the results are compared with those obtained from MD simulations using explicit water (EW) molecules in an all-atom representation. The initial structure of Dynorphin was assumed to be an alpha-helix between residues 1 and 9 as suggested from NMR measurements in micelles. The results obtained in the MD simulations show that the helical structure collapses early in the simulation, a behavior observed in the EW simulation and consistent with spectroscopic data that suggest that the peptide may adopt mainly an extended conformation in water. The dynamics of protein G calculated with the SCP implicit solvent model (SCP-ISM) reveals a stable structure that conserves all the elements of secondary structure throughout the entire simulation time. The average structures calculated from the trajectories with the implicit and explicit solvent models had a cRMSD of 1.1 A, whereas each average structure had a cRMSD of about 0.8A with respect to the X-ray structure. The main conformational differences of the average structures with respect to the crystal structure occur in the loop involving residues 8-14. Despite the overall similarity of the simulated dynamics with EW and SCP models, fluctuations of side-chains are larger when the implicit solvent is used, especially in solvent exposed side-chains. The MD simulation of Dynorphin was extended to 40 ns to study its behavior in an aqueous environment. This long simulation showed that the peptide has a tendency to form an alpha-helical structure in water, but the stabilization free energy is too weak, resulting in frequent interconversions between random and helical conformations during the simulation time. The results reported here suggest that the SCP implicit solvent model is adequate to describe electrostatic effects in MD simulation of both peptides and proteins using the same set of parameters. It is suggested that the present approach could form the basis for the development of a reliable and general continuum approach for use in molecular biology, and directions are outlined for attaining this long-term goal.

Comparative binding energy (COMBINE) analysis of OppA-peptide complexes to relate structure to binding thermodynamics

The periplasmic oligopeptide binding component (OppA) of the oligopeptide permease found in Gram-negative bacteria acts as a receptor for peptide transport across the cell membrane and is a potential target for antibacterial drug design. OppA exhibits broad specificity, binding to diverse peptides of 2-5 amino acid residues length. Crystallographic and calorimetric measurements have been carried out by Tame et al. of the binding of 28 peptides of sequence K-X-K to OppA, where X is a natural or nonnatural amino acid. Despite this extensive experimental characterization, a clear relationship between structural and thermodynamic parameters could not be readily identified, with a complicating factor being the observation of varying numbers of water molecules at the binding interface in the different complexes. Consequently, we have applied COMparative BINding Energy (COMBINE) analysis to derive quantitative structure-activity relationships (QSARs) for these 28 OppA-tripeptide complexes. This is the first application of COMBINE analysis to predict binding enthalpies and entropies, and predictive QSAR models were obtained for these quantities as well as for binding free energies. These QSAR models highlight several protein residues and bound water molecules in the binding site, as well as the electrostatic desolvation energies of the protein and the peptides, as responsible for most of the differences in binding thermodynamics between the peptides studied. The QSAR models aid rationalization of the determinants of binding affinity of the OppA:peptide complexes and provide guides for further ligand design. This study also points to the general applicability of COMBINE analysis to estimating thermodynamic parameters for protein-peptide complexes.

Efficient electrostatic solvation model for protein-fragment docking

A method is presented for the fast evaluation of the binding energy of a protein-small molecule complex with electrostatic solvation. It makes use of a fast preprocessing step based on the assumption that the main contribution to electrostatic desolvation upon ligand binding originates from the displacement of the first shell of water molecules. For a rigid protein, the precomputation of the energy contributions on a set of grids allows the estimation of the energy in solution of about 300 protein-fragment binding modes per second on a personal computer. The docking procedure is applied to five rigid binding sites whose size ranges from 17 residues to a whole protein of 107 amino acids. Using a library of 70 mainly rigid molecules, known micromolar inhibitors or close analogs are docked and prioritized correctly. The docking based rank-ordering of the library requires about 5 h and is proposed as a complementary approach to structure-activity relationships by nuclear magnetic resonance. Proteins 2001;42:256-268.

Calculating the local solvent chemical potential in crystal hydrates

Determining solvation patterns in biological systems is crucial in investigating the functional role water may play in structural stabilization and molecular recognition. Determining whether a particular position would be occupied by a solvent molecule requires the local thermodynamics to be known. In this work we introduce a simple and inexpensive approach based on grand canonical molecular simulations to determine the occupancy factors of the cavities. The method is applied to the test case of the sodium salt of hyaluronic acid. The results agree very well with experimental results and demonstrate the success of the method.

Theory of biomolecular recognition

Specific, noncovalent binding of biomolecules can only be understood by considering structural, thermodynamic, and kinetic issues. The theoretical foundations for such analyses have been clarified in the past year. Computational techniques for both particle-based and continuum models continue to improve and to yield useful insights into an ever wider range of biomolecular systems.

Comparative binding energy analysis of HIV-1 protease inhibitors: incorporation of solvent effects and validation as a powerful tool in receptor-based drug design

A comparative binding energy (COMBINE) analysis (Ortiz et al. J. Med. Chem. 1995, 38, 2681-2691) has been performed on a training set of 33 HIV-1 protease inhibitors, and the resulting regression models have been validated using an additional external set of 16 inhibitors. This data set was originally reported by Holloway et al. (J. Med. Chem. 1995, 38, 305-317), who showed the usefulness of molecular mechanics interaction energies for predicting the activity of novel HIV-1 protease inhibitors within the framework of the MM2X force field and linear regression techniques. We first used the AMBER force field on the same set of three-dimensional structures to check up on any possible force-field dependencies. In agreement with the previous findings, the calculated raw ligand-receptor interaction energies were highly correlated with the inhibitory activities (r2 = 0.81), and the linear regression model relating both magnitudes had an acceptable predictive ability both in internal validation tests (q2 = 0.79, SDEPcv = 0.61) and when applied to the external set of 16 different inhibitors (SDEPex = 1.08). When the interaction energies were further analyzed using the COMBINE formalism, the resulting PLS model showed improved fitting properties (r2 = 0.89) and provided better estimations for the activity of the compounds in the external data set (SDEPex = 0.83). Computation of the electrostatic part of the ligand-receptor interactions by numerically solving the Poisson-Boltzmann equation did not improve the quality of the linear regression model. On the contrary, incorporation of the solvent-screened residue-based electrostatic interactions and two additional descriptors representing the electrostatic energy contributions to the partial desolvation of both the ligands and the receptor resulted in a COMBINE model that achieved a remarkable predictive ability, as assessed by both internal (q2 = 0.73, SDEPcv = 0.69) and external validation tests (SDEPex = 0.59). Finally, when all the inhibitors studied were merged into a single expanded set, a new model was obtained that explained 91% of the variance in biological activity (r2 = 0.91), with very high predictive ability (q2 = 0.81, SDEPcv = 0.66). In addition, the COMBINE analysis provided valuable information about the relative importance of the contributions to the activity of individual residues that can be fruitfully used to design better inhibitors. All in all, COMBINE analysis is validated as a powerful methodology for predicting binding affinities and pharmacological activities of congeneric ligands that bind to a common receptor.

Solvation studies of DMP323 and A76928 bound to HIV protease: analysis of water sites using grand canonical Monte Carlo simulations

We examine the water solvation of the complex of the inhibitors DMP323 and A76928 bound to HIV-1 protease through grand canonical Monte Carlo simulations, and demonstrate the ability of this method to reproduce crystal waters and effectively predict water positions not seen in the DMP323 or A76928 structures. The simulation method is useful for identifying structurally important waters that may not be resolved in the crystal structures. It can also be used to identify water positions around a putative drug candidate docked into a binding pocket. Knowledge of these water positions may be useful in designing drugs to utilize them as bridging groups or displace them in the binding pocket. In addition, the method should be useful in finding water sites in homology models of enzymes for which crystal structures are unavailable.

Assessment of solvation effects on calculated binding affinity differences: trypsin inhibition by flavonoids as a model system for congeneric series

On the basis of molecular models of the interaction between trypsin and a series of seven structurally congeneric bioflavonoid inhibitors, the influence of solvation effects in the calculation of binding free energy differences in congeneric series has been assessed. The models were derived by making use of the X-ray crystal structure of bovine trypsin and the DOCK program, and the complementarity of the interactions between the functional groups of the docked molecules and the binding site region was corroborated independently with the GRID program. Interaction energies calculated for the complexes using molecular mechanics were found to correlate with the experimental inhibitory activities, although the quality of the correlation was dependent on the molecular modeling protocol. To understand why such correlations could be obtained in the absence of an explicit description of solvent effects, the in vitro activities were transformed into binding free energies, and continuum electrostatic theory was used to incorporate solvent effects by approximating them to the electrostatic contribution to the binding free energies. The results of our calculations show that, within this congeneric series, the cost in electrostatic free energy of desolvating both the enzyme binding site and the buried part of the inhibitors (delta Gdesolv) is roughly constant within the series. On the other hand, the electrostatic interaction energy (EeleLR) varies only slightly along the series in comparison with the van der Waals interaction (EVDWLR), and this variation is mostly solvent-independent, i.e., the reaction field energy of the solvent in the bound state (EsrfLR) makes almost a negligible contribution to the binding free energy differences. As a result, differences in binding free energy are dominated by the van der Waals term, while the electrostatic contribution is, to a good approximation, solvent-independent. A similar scenario may account for the good correlations frequently found between ligand activities and ligand-receptor interaction energies derived using plain molecular mechanics, although generality remains to be determined.

Analysis of the stability of looped-out and stacked-in conformations of an adenine bulge in DNA using a continuum model for solvent and ions

A combination of conformational search, energy minimization, and energetic evaluation using a continuum solvent treatment has been employed to study the stability of various conformations of the DNA fragment d(CGCAGAA)/d(TTCGCG) containing a single adenine bulge. The extra-helical (looped-out) bulge conformation derived from a published x-ray structure and intra-helical (stacked bulge base) model structures partially based on nuclear magnetic resonance (NMR) data were used as start structures for the conformational search. Solvent-dependent contributions to the stability of the conformations were calculated from the solvent exposed molecular surface area and by using the finite difference Poisson-Boltzmann approach. Three classes (I-III) of bulge conformations with calculated low energies can be distinguished. The lowest-energy conformations were found in class I, corresponding to structures with the bulge base stacked between flanking helices, and class II, composed of structures forming a triplet of the bulge base and a flanking base pair. All extra-helical bulge structures, forming class III, were found to be less stable compared with the lowest energy structures of class I and II. The results are consistent with NMR data on an adenine bulge in the same sequence context indicating an intra-helical or triplet bulge conformation in solution. Although the total energies and total electrostatic energies of the low-energy conformations show only relatively modest variations, the energetic contributions to the stability were found to vary significantly among the classes of bulge structures. All intra-helical bulge structures are stabilized by a more favorable Coulomb charge-charge interaction but destabilized by a larger electrostatic reaction field contribution compared with all extra-helical and most triplet bulge structures. Van der Waals packing interactions and nonpolar surface-area-dependent contributions appear to favor triplet class II structures and to a lesser degree also the intra-helical stacked bulge conformations. The large conformational variation found for class III conformers might add a favorable entropic contribution to the stability of the extra-helical bulge form.