Different microscopic and semimicroscopic approaches for calculations of electrostatic energies in macromolecules are examined. This includes the Protein Dipoles Langevin Dipoles (PDLD) method, the semimicroscopic PDLD (PDLD/S) method, and a free energy perturbation (FEP) method. The incorporation of these approaches in the POLARIS and ENZYME modules of the MOLARIS package is described in detail. The PDLD electrostatic calculations are augmented by estimates of the relevant hydrophobic and steric contributions, as well as the effects of the ionic strength and external pH. Determination of the hydrophobic energy involves an approach that considers the modification of the effective surface area of the solute by local field effects. The steric contributions are analyzed in terms of the corresponding reorganization energies. Ionic strength effects are studied by modeling the ionic environment around the given system using a grid of residual charges and evaluating the relevant interaction using Coulomb's law with the dielectric constant of water. The performance of the FEP calculations is significantly enhanced by using special boundary conditions and evaluating the long-range electrostatic contributions using the Local Reaction Field (LRF) model. A diverse set of electrostatic effects are examined, including the solvation energies of charges in proteins and solutions, energetics of ion pairs in proteins and solutions, interaction between surface charges in proteins, and effect of ionic strength on such interactions, as well as electrostatic contributions to binding and catalysis in solvated proteins. Encouraging results are obtained by the microscopic and semimicroscopic approaches and the problems associated with some macroscopic models are illustrated. The PDLD and PDLD/S methods appear to be much faster than the FEP approach and still give reasonable results. In particular, the speed and simplicity of the PDLD/S method make it an effective strategy for calculations of electrostatic free energies in interactive docking studies. Nevertheless, comparing the results of the three approaches can provide a useful estimate of the accuracy of the calculated energies. 0 1993 by John Wiley & Sons, Inc.

One of the most direct benchmarks for electrostatic models of macromolecules is provided by the pK a's of ionizable groups in proteins. Obtaining accurate results for such a benchmark presents a major challenge. Microscopic models involve very large opposing contributions and suffer from convergence problems. Continuum models that consider the protein permanent dipoles as a part of the dielectric constant cannot reproduce the correct self-energy. Continuum models that treat the local environment in a semi-microscopic way do not take into account consistently the protein relaxation during the charging process. This work describes calculations of pK a's in protein in an accurate yet consistent way, using the semi-microscopic version of the protein dipoles Langevin dipoles (PDLD) model, which treats the protein relaxation in the microscopic framework of the linear response approximation. This approach allows one to take into account the protein structural reorganization during formation of charges, thus reducing the problems with the use of the so-called “protein dielectric constant”, εp. The model is used in calculations of pK a's of the acidic groups of lysozyme, and the calculated results are compared to the corresponding results of discretized continuum (DC) studies. It is found that the present approach is more consistent than current DC models and also provides improved accuracies. Significant emphasis is given to the self-energy term, which has been pointed out in our early works but has been sometimes overlooked or presented as a small effect. The meaning of the dielectric constant εp used in DC models is clarified and illustrated, establishing the finding (e.g. King et. al., J. Phys. Chem. 1991, 95, 4366) that this parameter represents the contributions that are not treated explicitly in the given model, rather than the “true” dielectric constant. It is pointed out that recent suggestions to use large εp to obtain improved DC results might not be much different than our earlier suggestion to use a large effective dielectric for charge−charge interactions. This εp reduces the overestimate of charge−charge interactions relative to models that use small εp while not considering the protein relaxation explicitly. Unfortunately, the use of large εp does not reproduce consistently the self-energies of isolated ionized groups in protein interiors. The recent interest in taking protein flexibility into account in pK a calculations is addressed. It is pointed out that running MD over protein configurations will not by itself lead to a more consistent value of εp. It is clarified that a smaller value of εp, which is not really more (or less) consistent with the physics of the proteins, will be obtained if one uses our LRA (linear response approximation) formulation, generating configurations of both neutral and ionized states of the protein. It is also stated that such studies have been a standard part of our approach for some time. The present model involves a consecutive running of all-atom MD simulations of solvat...

The solution structures of the reduced and oxidized forms of the cytochrome c are used to reevaluate the reorganization energy for oxidation of cytochrome c. This is achieved by using the linear response approximation in concert with the NMR structures as pseudo energy constraints. Alternative estimates, obtained using a free energy perturbation approach employing umbrella sampling and a continuum dielectric approach, are also provided. The reorganization energy obtained is larger than that previously estimated using crystal structures of the protein. Nevertheless, the present estimate remains significantly smaller than the corresponding reorganization energy in water (9−15 kcal mol-1 as compared to ≈37 kcal mol-1 in water) and the protein contribution to the reorganization energy is only 8−10 kcal mol-1. This provides further support for the proposal that proteins assist in electron transfer reactions by reducing the relevant reorganization energies. The solution structures are also used to estimate the redox potential of cytochrome c. Several strategies are employed including a newly formulated scaled linear response approximation. The calculations agree reasonably well with the observed redox potential. Analysis of the group contributions to the reorganization energy and redox potential reveals a clear energetic linkage between these fundamental parameters of electron transfer and a redox-dependent surface feature likely to influence recognition of cytochrome c by its redox partners. Specifically, the rearrangement of Ile81 and other residues at the heme edge upon a change in oxidation state gives rise to a large contribution to both the redox potential and the reorganization energy. Finally this work is used to explore and illustrate the meaning of macroscopic dielectric models. It is shown that the “proper” dielectric constant depends strongly on the model used since it basically represents the implicit contributions of the given model rather than a fundamental physics. Thus we obtain different effective dielectric constants for different treatments of redox potential and reorganization energy.

Several strategies for evaluation of the protein-ligand binding free energies are examined. Particular emphasis is placed on the Linear Response Approximation (LRA) (Lee et. al., Prot Eng 1992;5:215-228) and the Linear Interaction Energy (LIE) method (Aqvist et. al., Prot Eng 1994;7:385-391). The performance of the Protein Dipoles Langevin Dipoles (PDLD) method and its semi-microscopic version (the PDLD/S method) is also considered. The examination is done by using these methods in the evaluating of the binding free energies of neutral C2-symmetric cyclic urea-based molecules to Human Immunodeficiency Virus (HIV) protease. Our starting point is the introduction of a thermodynamic cycle that decomposes the total binding free energy to electrostatic and non-electrostatic contributions. This cycle is closely related to the cycle introduced in our original LRA study (Lee et. al., Prot Eng 1992;5:215-228). The electrostatic contribution is evaluated within the LRA formulation by averaging the protein-ligand (and/or solvent-ligand) electrostatic energy over trajectories that are propagated on the potentials of both the polar and non-polar (where all residual charges are set to zero) states of the ligand. This average involves a scaling factor of 0.5 for the contributions from each state and this factor is being used in both the LRA and LIE methods. The difference is, however, that the LIE method neglects the contribution from trajectories over the potential of the non-polar state. This approximation is entirely valid in studies of ligands in water but not necessarily in active sites of proteins. It is found in the present case that the contribution from the non-polar states to the protein-ligand binding energy is rather small. Nevertheless, it is clearly expected that this term is not negligible in cases where the protein provides preorganized environment to stabilize the residual charges of the ligand. This contribution can be particularly important in cases of charged ligands. The analysis of the non-electrostatic term is much more complex. It is concluded that within the LRA method one has to complete the relevant thermodynamic cycle by evaluating the binding free energy of the "non-polar" ligand, l;, where all the residual charges are set to zero. It is shown that the LIE term, which involves the scaling of the van der Waals interaction by a constant beta (usually in the order of 0.15 to 0.25), corresponds to this part of the cycle. In order to elucidate the nature of this non-electrostatic term and the origin of the scaling constant beta, it is important to evaluate explicitly the different contributions to the binding energy of the non-polar ligand, DeltaG(bind,l;). Since this cannot be done at present (for relatively large ligands) by rigorous free energy perturbation approaches, we evaluate DeltaG(bind,l;) by the PDLD approach, augmented by microscopic calculations of the change in configurational entropy upon binding. This evaluation takes into account the van der Waals, hydrophobic, water penetration and entr...

No abstract

The quantum dynamics of the primary photoisomerization event in bacteriorhodopsin is studied by a semiclassical trajectory approach. The relevant surface crossing probability is evaluated from the wave functions and potential surfaces of a hybrid quantum mechanical/molecular mechanics (QM/MM) Hamiltonian of the complete chromophore-protein-solvent system. The QM/MM model combines consistently the quantum mechanical Hamiltonian of the chromophore with the microscopic electric field of the ionized groups and induced dipoles of the protein-solvent system. The QCFF/PI Hamiltonian of the chromophore is adjusted to reproduce relevant ab initio results. The nonadiabatic coupling term <ψ 1 |∂ψ 0 /∂t> is calculated numerically from the corresponding wave functions. The simulations are performed by combining the ENZYMIX and QCFF/PI molecular modeling programs. The effect of the protein on the absorption spectrum of the chromophore is examined. It is found that this spectrum reflects the effect of the protein permanent dipoles, ionized residues, water molecules (in and around the protein), and the induced dipoles of the protein plus water system. Next, we probe the motion along the excited state surface. It is demonstrated, in agreement with our early study and more recent works, that the motion starts with bond vibrations and evolves to a torsional motion. It is also found that we are dealing with an overdumped motion. Major emphasis is placed on the nature of the surface crossing process. In particular, we try to examine the origin of the very large probability of crossing in the π/2 region. A large crossing probability was obtained first in our early simulation (Warshel, A. Nature 1976, 260, 679), but its origin was not explored in details. Such large crossing probabilities can be obtained by passing through strict conical intersections (where the two surfaces "touch" each other) or by passing through regions with large nonadiabatic coupling and small energy gap (such regions are usually close to conical intersections). It is found that some trajectories pass through strict conical intersections, whereas others cross through regions with nonzero energy gap and a large nonadiabatic coupling. This feature helps probably to ensure the stability of the photobiological process with regards to various mutations. The average surface crossing probability and our previously derived expression (Weiss, R. M.; Warshel, A. J. Am. Chem. Soc. 1979, 101, 6131) appear to provide an excellent approximation for the calculated quantum yield. Furthermore, the calculated quantum yield reproduces the corresponding observed value. Finally, we examine the behavior of trajectories that cross to the ground state before the π/2 region. Our finding that these trajectories are deflected backward allow us to exclude models where the surface crossing occurs before the π/2 region.

One of the fundamental challenges in biotechnology and in biochemistry is the ability to design effective enzymes. Doing so would be a convincing manifestation of a full understanding of the origin of enzyme catalysis. Despite an impressive progress, most of the advances on this front have been made by placing the reacting fragments in the proper places, rather than by optimizing the environment preorganization, which is the key factor in enzyme catalysis. Rational improvement of the preorganization would require approaches capable of evaluating reliably the actual catalytic effect. This work takes apreviously designed kemp eliminases as a benchmark for a computer aided enzyme design, using the empirical valence bond as the main screening tool. The observed absolute catalytic effect and the effect of directed evolution are reproduced and analyzed (assuming that the substrate is in the designed site). It is found that, in the case of kemp eliminases, the transition state charge distribution makes it hard to exploit the active site polarity, even with the ability to quantify the effect of different mutations. Unexpectedly, it is found that the directed evolution mutants lead to the reduction of solvation of the reactant state by water molecules rather that to the more common mode of transition state stabilization used by naturally evolved enzymes. Finally it is pointed out that our difficulties in improving Kemp eliminase are not due to overlooking exotic effect, but to the challenge in designing a preorganized environment that would exploit the small change it charge distribution during the formation of the transition state.computer aided enzyme design | empirical valence bond | directed evolution R ational enzyme design is expected to have a great potential in industrial application and eventually in medicine (1). Furthermore, the ability to design efficient enzymes might be considered as the best manifestation of a true understanding of enzyme catalysis. However, at present there has been a limited success in most attempts of rational enzyme design, and the resulting constructs have been much less effective than the corresponding natural enzymes (1). Furthermore, despite the progress in directed evolution (e.g., ref.2), we do not have unique rationales for the resulting rate enhancements.Most attempts to identify the problems with the current rational design approaches (for review, see ref. 1) have not been based on actual simulations of the given effect. In fact, it has been argued (3,4), that the problems are due to the incomplete modeling of the transition state (TS) and to the limited awareness to the key role of the reorganization energy. Even a recent attempt to use a molecular orbital-combined quantum mechanical /molecular mechanics (MO-QM/MM) approach (5) has not provided a reasonable estimate of the observed catalytic effect or the trend of the mutational effects in an artificially design enzyme. Thus, reproducing the effect of directed evolution and eventually obtaining better performance in enzyme de...

The idea that enzymes accelerate their reactions by entropic effects has played a major role in many prominent proposals about the origin of enzyme catalysis. This idea implies that the binding to an enzyme active site freezes the motion of the reacting fragments and eliminates their entropic contributions, (⌬S cat ‡ )Ј, to the activation energy. It is also implied that the binding entropy is equal to the activation entropy, (⌬S w ‡ )Ј, of the corresponding solution reaction. It is, however, difficult to examine this idea by experimental approaches. The present paper defines the entropic proposal in a rigorous way and develops a computer simulation approach that determines (⌬S ‡ )Ј. This approach allows us to evaluate the differences between (⌬S ‡ )Ј of an enzymatic reaction and of the corresponding reference reaction in solution. Our approach is used in a study of the entropic contribution to the catalytic reaction of subtilisin. It is found that this contribution is much smaller than previously thought. This result is due to the following: (i ) Many of the motions that are free in the reactants state of the reference solution reaction are also free at the transition state. Many prominent proposals (e.g., see refs. 1 and 2) and textbooks that consider biochemical systems (e.g., refs. 3 and 4) invoke entropic contributions as major factors in enzyme catalysis. These proposals, which are intuitively very appealing (e.g., see ref. 5), have assumed that the large configurational space available for the reacting fragments in water would be drastically restricted in the enzyme active site. It has been thus deduced that this should lead to large entropic contributions to the difference between the activation barrier in the enzyme and in the reference solution reaction. However, the validity of these proposals is far from being obvious (6, 7). For example, the very inf luential proposal introduced by Page and Jencks (1) ref lects the assumption that the formation of the transition state in a bimolecular reaction in solution involves complete loss of three translational and three rotational degrees of freedom. However, two or more of these degrees of freedom are usually almost free in the transition state (see below). More serious is the implicit assumption that the entropic contribution to catalysis is given approximately by the negative of the binding entropy (see below). Other problems with simple estimates of the entropic contribution will be mentioned in the next section.The main stumbling block for determining the validity of the entropic proposal is the absence of direct experimental information about the corresponding contribution of the reacting fragments to the activation entropy in the enzyme and in solution. In this respect it is interesting to note the recent analysis of cytidine deaminase by Wolfenden and co-workers (8). This study found that the entropies of activation in the enzyme and in water are very similar and that the overall catalysis is due to enthalpic effects. Interestingly, it was found that ...