3168
J. Phys. Chem. B 2008, 112, 3168-3176
Hepatitis C Virus NS5B Polymerase: QM/MM Calculations Show the Important Role of the Internal Energy in Ligand Binding Jerry M. Parks,† Rama K. Kondru,‡ Hao Hu,† David N. Beratan,† and Weitao Yang*,† Department of Chemistry, Duke UniVersity, 124 Science DriVe, 5301 French Science Center, Durham, North Carolina 27708-0346, and Roche Palo Alto LLC, 3431 HillView AVenue, Palo Alto, California 94304 ReceiVed: August 28, 2007; In Final Form: December 12, 2007
The inter- and intramolecular interactions that determine the experimentally observed binding mode of the ligand (2Z)-2-(benzoylamino)-3-[4-(2-bromophenoxy)phenyl]-2-propenoate in complex with hepatitis C virus NS5B polymerase have been studied using QM/MM calculations. DFT-based QM/MM optimizations were performed on a number of ligand conformers in the protein-ligand complex. Using these initial poses, our aim is 2-fold. First, we identify the minimum energy pose. Second, we dissect the energetic contributions to this pose using QM/MM methods. The study reveals the critical importance of internal energy for the proper energy ranking of the docked poses. Using this protocol, we successfully identified three poses that have low RMSD with respect to the crystallographic structure from among the top 20 initially docked poses. We show that the most important energetic component contributing to binding for this particular protein-ligand system is the conformational (i.e., QM internal) energy.
I. Introduction Understanding the intramolecular and intermolecular effects that govern protein-ligand binding is an important goal of rational drug design. Molecular recognition involves a complex interplay between numerous covalent and noncovalent interactions. Ligand binding may be influenced by hydrogen bonding, electrostatic interactions, steric effects, solvation/desolvation, cation-π interactions, conformational energetics, or a combination of these effects. The relative magnitudes of individual effects vary from protein to protein and from ligand to ligand. Thus, it is often necessary to examine energetic contributions to binding in protein-ligand systems on a case-by-case basis. Energy decomposition analysis affords the possibility of determining which individual energetic components are most important to the binding process. Over the past three decades, quantum mechanics/molecular mechanics (QM/MM) methods have been used to solve numerous problems in the simulation of biological systems.1 The basic premise of the QM/MM approach is that a small portion of a macromolecular system is described quantum mechanically while the remainder of the system is modeled by molecular mechanics energy functions.2 Our lab has developed ab initio and density functional theory (DFT) QM/MM methods primarily for the study of enzymatic reaction mechanisms.3-10 One limitation of standard, nonpolarizable MM-based methods is the neglect of electronic polarization effects exerted on a ligand by the protein environment and vice versa. These charge effects can be important for correctly describing electrostatic complementarity in the binding pocket. QM calculations include polarization effects explicitly by adding an electrostatic term to the Hamiltonian that represents the protein environment as a set of point charges. * To whom correspondence should be addressed. † Duke University. ‡ Roche Palo Alto LLC.
Quantum chemical calculations are capable of providing accurate molecular geometries and conformational energetics, and are able to describe ligand- (and to a limited extent, protein-) charge polarization, but they have greater computational demands than purely classical approaches. DFT methods are frequently used in QM/MM calculations involving enzymatic systems. Commonly used functionals, including the B3LYP exchange-correlation functional,11,12 generally provide a reasonable balance between computational efficiency and accurate molecular geometries and energetics, although certain pathologies do exist.13-15 Among the well-known deficiencies of DFT methods is their poor treatment of dispersion interactions.16,17 DFT functionals lack the inclusion of instantaneously fluctuating dipole interactions, and thus do not include van der Waals effects explicitly. However, it is quite common in QM/MM calculations to describe the van der Waals interactions between the QM subsystem and nearby MM atoms via a nonbonded LennardJones term.1 Wu et al.18 devised an empirical van der Waals correction for use with standard DFT functionals. An attractive R-6 term was added to the total energy, yielding results that were comparable with MP2 calculations. Jurecˇka et al.19 used a similar concept and augmented various DFT functionals with a damped empirical dispersion term to achieve good results on a test set of 58 van der Waals complexes. The choice of VDW parameters is somewhat arbitrary because atom types must be defined for QM atoms by analogy. However, because of the short-range nature of VDW interactions, the exact values of the parameters are not critical.20,21 Thus, the standard approach of using MM VDW parameters for the description of QM/MM VDW interactions is at least somewhat justified and has been shown to be acceptable in practical calculations. Hepatitis C Virus Polymerase. Hepatitis C affects millions of people worldwide and often results in chronic infection and
10.1021/jp076885j CCC: $40.75 © 2008 American Chemical Society Published on Web 02/14/2008
Hepatitis C Virus NS5B Polymerase
Figure 1. (2Z)-2-(benzoylamino)-3-[4-(2-bromophenoxy)phenyl]-2propenoate ligand.
J. Phys. Chem. B, Vol. 112, No. 10, 2008 3169 (2Z)-2-(benzoylamino)-3-[4-(2-bromophenoxy)phenyl]-2-propenoate (Figure 1). The crystallographic structure of this protein-ligand complex is known, so we first validate our QM/MM methodology by optimizing the ligand as found in the crystal structure and comparing it with the original coordinates. Next, we generate several initial docked poses with the program FlexX.23-25 Using these initial poses, we first identify the minimum energy pose. Next, we dissect the energetic contributions to these poses using QM/MM methods. Thus, the goal of this study is to answer two questions: Can the QM/MM energy alone provide a meaningful scoring function for the HCV protein-ligand system and reproduce the crystallographic binding mode? And if so, why does this particular ligand bind in its crystallographic pose? Theory. QM/MM Energy Expression. For our QM/MM calculations, we adopt the same formalism used by Zhang et al.,3 although the pseudobond method is not employed because there are no covalent bonds between the protein and ligand in the HCV system. If necessary, protein residues could be easily included in the QM subsystem using the pseudobond method, although the computational cost would increase significantly. In our approach, the QM/MM total energy expression is partitioned as
Etotal(rQM, rMM) ) EQM(rQM) + EQM/MM(rQM, rMM) + EMM(rMM) (1)
Figure 2. Ligand and adjacent binding site residues of HCV polymerase. The crystallographic water molecule in the binding site is also shown. (Figure generated using VMD.43,44)
chronic liver disease. There is no vaccine for hepatitis C, and combination therapies including interferon and ribavirin are not entirely successful in treating the disease. The HCV NS5B RNA-dependent RNA polymerase enzyme is required for viral replication. As a result, this enzyme represents a promising potential target for small-molecule inhibitors. Pfefferkorn et al.22 obtained a 2.5 Å crystal structure of HCV NS5B polymerase with the ligand molecule shown in Figure 1. They observed several notable interactions between the protein and ligand, including a water-mediated hydrogen bond between the carboxylate group of the ligand and the backbone N-H of Ser556. A hydrogen bond is present between the carbonyl of the benzamide group and the backbone N-H of Tyr448. They also found that the A and B rings of the ligand are located in a hydrophobic region consisting of Pro197, Leu384, Met414, Tyr415, and Tyr448. The ligand and surrounding residues in the binding site are displayed in Figure 2. Goals of This Study. In this work, we have applied a straightforward QM/MM approach as a first step in understanding ligand binding in HCV NS5b polymerase. A more robust protocol, including a more realistic description of the solvent, a more accurate description of the electrostatic interactions, protein side chain flexibility, and dynamic sampling would undoubtedly provide a more detailed description of specific protein-ligand interactions. However the computational cost of such a simulation would be significantly higher. We aim here to use a standard scoring function-based docking program to generate initial ligand conformers and then optimize the resulting protein-ligand complexes using QM/MM calculations. The system we have chosen is the C∆21 subunit of hepatitis C virus NS5B polymerase complexed with the ligand
where rQM and rMM represent the coordinates of the QM and MM subsystems, respectively. The second term on the RHS, EQM/MM, refers to the interaction energy between the two subsystems. EQM/MM is decomposed to yield
EQM/MM(rQM, rMM) ) Eelec QM/MM(rQM, rMM) + cov Evdw QM/MM(rQM, rMM) + EQM/MM(rQM, rMM) (2) elec where EQM/MM is the electrostatic interaction energy between vdw is the van der Waals the QM and MM subsystems, EQM/MM cov interaction energy between the two subsystems, and EQM/MM includes any covalent bonding interactions between the two subsystems. Often, as in the present case, there are no covalent linkages between the protein and the ligand, so the last term is zero. elec In practice, we combine EQM and EQM/MM because it is necessary to compute these terms simultaneously using a quantum chemistry software package. This effectiVe Hamiltonian term is obtained by performing a QM calculation in which the molecule is surrounded by a set of external point charges. The point charges are centered at the nuclear positions of the MM atoms within a predefined cutoff distance. The magnitude of each charge is simply the partial charge from the MM force field. Thus, the incorporation of electrostatic polarization of the QM subsystem by the surrounding protein environment is accomplished in a direct way. For the purpose of decomposing the energy to obtain the electrostatic interactions between the QM and MM subsystems, the QM subsystem is approximated by a set of charges obtained via electrostatic potential (ESP) fitting. ESP charges are computed by fitting partial charges centered on the QM nuclei so that they reproduce the molecular electrostatic potential (MEP). In this way, the classical Coulombic interaction energy between the QM and MM subsystems can be evaluated rapidly. The relationship between the energy of the effective Hamiltonian and its two components, namely the polarized internal energy and the classical QM/MM interaction energy is given by
3170 J. Phys. Chem. B, Vol. 112, No. 10, 2008 elec 〈Ψ|Heff|Ψ〉 ) EQM + Eelec QM/MM ) E′QM + E′QM/MM
Parks et al.
(3)
where elec ) E′QM/MM
∑ ∑ i∈QM j∈MM
Qiqj rij
(4)
Here Qi is the ESP charge on QM atom i and qj is the MM point charge on atom j. We may then rearrange eq 3 to yield elec E′QM ) (EQM + Eelec QM/MM) - E′QM/MM
(5)
where E′QM is the QM energy in the presence of the MM point charges with the classical Coulombic interaction energy subtracted. Thus, E′ is the internal energy of the QM subsystem that is polarized by the MM environment. To model the effects of solvation, we introduce a distancedependent dielectric constant26-29 in the expression for the electrostatic interaction energy elec,(r) ) E′QM/MM
∑ ∑
i∈QM j∈MM
1 Qiqj (r) rij
(6)
Here, (r) is the dielectric constant, where we have chosen two simple forms, (r)) rij or (r) ) 4rij. Thus, we mimic bulk solvation in a simple way by screening the charge-charge interactions between the protein and ligand. Note that the unscreened MM point charges are used in the effective Hamiltonian during the QM/MM optimization, but they are not included in the QM/MM total energy expression. Instead, the charges screened by the dielectric are used to compute the final QM/MM energy. This was done for simplicity of implementation. It proved unnecessary in the present case to include screening in the optimization, although such a modification would be desirable and is likely to yield improved results. vdw The van der Waals contributions (EQM/MM ) to the nonbonded interaction energy between the two subsystems are treated in a purely classical manner using attractive and repulsive LennardJones terms. As in the MM subsystem, van der Waals parameters from the MM force field are assigned to each atom in the QM subsystem. Atom types are assigned to the QM atoms by analogy (from the MM force field) and are the only parameters required for the ligand molecule. Finally, EMM consists of the bond, angle, and torsional terms as in standard MM force fields. However, since the protein structure is identical in each of the protein-ligand complex structures, this term is ignored. As such, our expression for the QM/MM energy is elec,(r) vdw Etotal QM/MM ) E′QM + E′QM/MM + EQM/MM
(7)
II. Methods Initial Docking Procedure. The crystal structure of the C∆21 subunit of NS5B HCV polymerase was obtained from the Protein Data Bank PDB code 1YVF.22 The program FlexX23-25 was used for the preliminary docking step. The binding environment was defined to consist of all protein atoms in the crystal structure that were within 6.0 Å of the ligand based on the crystal structure position, and all crystallographic waters were removed. Partial charges were assigned by FlexX and 200 poses were generated. All FlexX parameters were unchanged from the default values. No attempt was made to refine the empirical parameters or charges because we simply used FlexX
to generate initial conformers for QM/MM optimization. The top 20 poses were retained for further optimization. System Preparation. Hydrogens were added to the protein using the H++ server.30-32 In this approach, the PoissonBoltzmann equation is solved within a continuum dielectric representation of the solvent to compute the pKa values of ionizable protein sidechains and assign protonation states. An internal dielectric of 6.0, external dielectric of 80.0, and ionic strength of 0.15 M were used. All histidines were predicted to be neutral, and hydrogens were placed at the position by the H++ server. Cys14 was predicted to be anionic and was modeled as such, although this residue is far from the binding pocket and does not contribute to binding. All other ionizable sidechains were modeled in their standard pH 7 protonation states. QM/MM Procedure. The 43 atoms of the ligand were chosen as the QM subsystem. All other atoms were represented by molecular mechanics using the Charmm22 force field.33 Nonbonded van der Waals parameters for Br were obtained from the General Amber Force Field (GAFF) set.34 QM/MM calculations were performed with the program Sigma9,35,36 interfaced with a modified version of Gaussian 03.37 The B3LYP exchange-correlation functional11,12 and the 6-31G* basis set were used for all QM calculations, and the Merz-SinghKollman38 scheme was used for ESP charge fitting. For the interactions between the QM and MM subsystems, a nonbonded cutoff of 15.0 Å was used. Loose convergence criteria were used for the optimizations of the QM subsystems. Initially, an MM optimization of the hydrogen atoms of the protein was performed using conjugate gradient minimization. From the MM docking results using FlexX, the top 20 conformers were each placed into the binding pocket of identically prepared protein structures. III. Results Initial Conformer Generation. In the following sections, the numbering of each conformer indicates its rank as predicted by FlexX. Thus, conformer 1 was the highest scoring pose according to the FlexX scoring function. Although the RMSD for conformer 1 relative to the crystal structure conformer was more than 7 Å, three of the 20 highest ranking poses predicted by FlexX were near-native. Conformers 3, 4 and 12 had RMSD values relative to the crystal structure of less than 1.0 Å. The 20 highest ranking ligand conformers from FlexX were then subjected to QM/MM optimization. We also performed the same optimization for the unmodified crystal structure complex, which was used as a reference for both energetic and structural comparisons. It is important to note that the protein structures for each conformer were identical. The protein hydrogens were optimized in a separate, purely MM optimization with no ligand present. This protein structure was then replicated and each ligand conformer was implanted in the binding pocket. The heavy atoms of the protein in both the FlexX and QM/MM optimizations were identical, so no fitting of the ligand structures was necessary. In each QM/MM optimization, the protein was held fixed. This procedure allowed a direct comparison of the energies of each complex. If the protein were allowed to move, it is likely that the different complexes would end up in irrelevant local minima. Exclusion of Water Molecules in the Binding Pocket. Another important point to address is the exclusion of solvent molecules from our simulations. The crystal structure22 reveals a water molecule involved in a hydrogen bond network between the carboxylate of the ligand and the backbone N-H of Ser556
Hepatitis C Virus NS5B Polymerase
J. Phys. Chem. B, Vol. 112, No. 10, 2008 3171
TABLE 1: QM/MM Energies and Individual Interaction Terms Using a Distance-Dependent Dielectric Constant of E(r) ) r or E(r) ) 4ra conformer
Egas
E′
′elec,(r))r EQM/MM
′elec,(r))4r EQM/MM
vdw EQM/MM
total,)r EQM/MM
total,)4r EQM/MM
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
9.18 9.19 -0.88 -0.87 19.47 19.47 16.01 16.01 10.82 10.86 10.37 1.72 17.08 17.09 12.22 11.92 16.88 15.99 15.99 15.77
23.86 23.85 0.16 0.18 38.28 38.28 30.96 30.96 26.27 26.31 24.83 1.57 31.36 31.35 25.00 24.56 31.09 29.96 29.94 30.61
-29.74 -29.72 -1.65 -1.66 -37.91 -37.91 -29.70 -29.71 -35.15 -35.21 -37.26 -0.30 -34.68 -34.63 -29.44 -28.83 -30.16 -28.23 -28.22 -29.36
-7.44 -7.43 -0.41 -0.41 -9.48 -9.48 -7.43 -7.43 -8.79 -8.80 -9.31 -0.07 -8.67 -8.66 -7.36 -7.21 -7.54 -7.06 -7.06 -7.34
21.40 21.39 10.75 10.70 26.79 26.80 23.18 23.18 22.69 22.72 22.86 5.01 18.77 18.66 20.17 20.04 24.44 22.79 22.78 22.30
15.52 15.52 9.27 9.23 27.16 27.17 24.43 24.43 13.81 13.82 10.44 6.29 15.45 15.37 15.72 15.77 25.37 24.52 24.50 23.56
37.82 37.81 10.50 10.47 55.59 55.60 46.71 46.72 40.17 40.23 38.38 6.51 41.46 41.35 37.80 37.40 47.99 45.69 45.67 45.58
a All values are relative to the QM/MM optimized crystal structure and are in kcal mol-1. For the crystal structure, the gas-phase energy of the ligand, Egas ) -2368710.27; the polarized QM internal energy (eq 5), E′ ) -2368707.74; the QM/MM electrostatic interaction energy (eq 6), ′elec,)r elec,)4r vdw total,)r EQM/MM ) -29.02, and EQM/MM ) -7.25; the van der Waals interaction energy, EQM/MM ) -43.78; and the QM/MM total energy (eq 7), EQM/MM total,)4r ) -2368780.53, and EQM/MM ) -2368758.77.
(see Figure 2). We chose to remove the water molecule because positions of water molecules in a binding pocket are not generally known a priori in docking studies. Inclusion of this water molecule would most likely improve our results because it would facilitate the hydrogen bond bridging interaction involving Ser556, although this assertion was not investigated. RMSD of Initial Versus QM/MM-Docked Conformers. To demonstrate how much each structure changed as a result of the QM/MM optimization procedure, we calculated the RMSD values between the 20 initial docking predictions from FlexX and their corresponding QM/MM optimized structures. The same procedure was performed for the crystal structure, except that the initial conformation of the crystal structure was simply based on the crystallographic coordinates. The heavy atom RMSD between the pre- and post-QM/MM optimized ligand from the crystal structure was only 0.21 Å. This value is quite encouraging considering the resolution of the crystal structure and that no explicit solvent molecules were present in the optimization. Thus, our exclusion of the active site water molecule appears to be acceptable in this case. This analysis also serves to validate the level of QM theory we have chosen. For all other conformers, it was found that the positions of the heavy atoms changed more significantly during the optimization. The RMSD values ranged from a relatively small 0.42 Å for conformers 7 and 8, to over 1.4 Å for conformer 11. These results indicate that the conformations assigned by FlexX did not correspond to minima on the potential energy surface of the system. This is to be expected considering the use of default charges and parameters in the generation of the initial conformers. QM/MM Energetics. As a means of approximating bulk solvation, we found the use of a distance-dependent dielectric constant to be essential. This simple approximation has been used successfully in many studies26-29 and serves to screen the bare charge-charge interactions that would otherwise be screened by solvent molecules. In this way, we reduced the errors resulting from artificially favorable interactions between the ligand and certain protein sidechains that would otherwise be strongly solvated in solution. It was found from molecular
dynamics simulations in explicit water that the region near the ligand carboxylate in the bound complex is quite solventaccessible and hydrophilic (data not shown). Note that the distance-dependent dielectric constant was not used in the optimization; it was used only in the final QM/MM energy evaluation. The QM/MM total energy of each complex was decomposed into individual components to provide a detailed comparison of the various conformers. In this way, we were also able to determine which energetic terms dominate in determining the experimental binding mode. The energies were evaluated with a distance-dependent dielectric constant of (r) ) rij and (r) ) 4rij. In each case, the energies were decomposed into contributions from the polarized QM internal energy, E′, the elec,(r) classical QM/MM electrostatic interaction energy, E′QM/MM , vdw and the van der Waals energy, EQM/MM. Both sets of energies are shown in Table 1. In both sets of calculations, the three lowest energy complexes corresponded to conformers 12, 3, and 4, thus providing a consistent picture between the two variants even though the overall magnitude of the relative energies is significantly different for the two dielectric models. From Table 1, we see that conformer 12 is 6.29 kcal mol-1 higher in energy than the crystal structure for (r) ) rij, and conformers 3 and 4 are around 9.2 kcal mol-1 higher than the crystallographic binding mode. Conformer 1, which was predicted to be the correct binding mode in the FlexX docking calculations, was 15.5 kcal mol-1 less stable than the crystal structure for ) rij. RMSD Versus the Crystal Structure after QM/MM Optimization. To determine whether the QM/MM energy of the docked complexes was sufficient to predict the correct docking conformation, we calculated the heavy atom RMSD of each ligand conformer relative to the ligand in the crystal structure. Conformers 3 and 4 were essentially identical and should be regarded as such. Conformers 12, 3, and 4 had RMSD values of 0.92, 0.87, and 0.87 Å, respectively. In addition to having the lowest QM/MM energy, conformer 12 had the desirable quality of an RMSD less than 1.0 Å. Conformers 3 and 4 were slighly more similar to the crystal structure in terms of RMSD values but were approximately 4 kcal mol-1 higher in energy
3172 J. Phys. Chem. B, Vol. 112, No. 10, 2008
Parks et al.
Figure 3. Comparison of energies of conformers 1-20 with their corresponding RMSD values. All energies are relative to the QM/MM optimized elec,)4r vdw crystal structure (in kcal mol-1). Top panel: QM/MM total energies and individual energetic terms. Red ) E′,cyan ) E′QM/MM , blue ) EQM/MM , total . Middle panel: RMSD of QM/MM optimized conformers relative to the crystal structure. Bottom panel: Gas-phase QM energies black )EQM/MM vdw at the optimized QM/MM geometry. Black ) Gas-phase single-point energies (Egas), blue ) Egas + EQM/MM (Color available online).
than conformer 12. The most obvious structural difference between conformer 12 and conformers 3 and 4 was that the phenyl substituent of the benzamide group was rotated slightly in conformer 12 relative to the other two poses (Figure 4). This difference is attributed predominantly to steric (i.e., van der Waals) effects, although the QM internal energy and polarization play a small role. All other complexes had significantly higher energies and the ligand poses differed considerably from the experimental binding mode with RMSD values between 6 and 8 Å. In fact, many of the conformers occupied a predominantly hydrophilic region adjacent to the binding pocket. A comparison of the QM/MM energies and the RMSD values for the 20 QM/MM optimized conformers is presented in Figure 3. The top panel illustrates the partitioning of the individual energetic terms from the QM/MM optimization, and the middle panel confirms the accuracy of the QM/MM energies by showing that the lowest energy structures are conformationally most similar to the crystal structure. Although there appears to be a correlation between the RMSD and QM/MM energy, it is only meaningful in that the lowest energy structures yield the smallest RMSD values. Thus, conformer 12 can be regarded as our best guess of the correct binding mode based purely on the QM/MM energy. Figure 4 shows conformer 12 superimposed on the crystallographic ligand conformer. It is evident that the overall orientation observed in the crystal structure has been reproduced by our optimization method. As previously noted, no explicit water molecules were included in the optimization, which is likely a dominant source of error in the calculations. We should also consider the possibility that conformers 3 and 4 represent another accessible binding conformation based on their similar but slightly higher energy. It is quite possible that the systematic errors in the QM/MM optimization procedure could affect the absolute ordering of the conformer energies. However, conformers 3, 4 and 12 occupy the same binding site and are quite similar in structure.
Figure 4. Comparison of (a) the predicted binding mode (conformer 12) and (b) conformers 3 and 4 from QM/MM optimization superimposed on the ligand in the crystal structure (dark gray). The heavyatom RMSD was 0.92 Å for conformer 12 and 0.87 Å for conformers 3 and 4.
Gas-Phase Energetics. Gas-phase single-point QM energies were calculated for each conformer to understand the conformational energetics of the ligand molecule. Interestingly, the overall trend in the ordering of the gas-phase energies roughly
Hepatitis C Virus NS5B Polymerase
J. Phys. Chem. B, Vol. 112, No. 10, 2008 3173
TABLE 2: Comparison of Classically Computed QM/MM Electrostatic Interaction Energies Using Gas-Phase QM ESP Charges and QM/MM ESP Chargesa conformation
Eelec,)r gas
elec,)r EQM/MM
Eelec,)r pol
Eelec,)4r gas
elec,)4r EQM/MM
Eelec,)4r pol
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
-18.22 -18.22 -1.10 -1.09 -24.22 -24.22 -20.30 -20.29 -23.08 -23.12 -27.79 -1.13 -23.72 -23.70 -21.54 -21.04 -21.99 -18.99 -18.99 -19.84
-29.74 -29.72 -1.65 -1.66 -37.91 -37.91 -29.70 -29.71 -35.15 -35.21 -37.26 -0.30 -34.68 -34.63 -29.44 -28.83 -30.16 -28.23 -28.22 -29.36
-11.53 -11.51 -0.56 -0.56 -13.69 -13.69 -9.41 -9.42 -12.07 -12.09 -9.47 0.83 -10.95 -10.93 -7.90 -7.79 -8.17 -9.24 -9.23 -9.51
-4.55 -4.55 -0.27 -0.27 -6.06 -6.05 -5.07 -5.07 -5.77 -5.78 -6.95 -0.28 -5.93 -5.93 -5.39 -5.26 -5.50 -4.75 -4.75 -4.96
-7.44 -7.43 -0.41 -0.41 -9.48 -9.48 -7.43 -7.43 -8.79 -8.80 -9.31 -0.07 -8.67 -8.66 -7.36 -7.21 -7.54 -7.06 -7.06 -7.34
-2.88 -2.88 -0.14 -0.14 -3.42 -3.42 -2.35 -2.35 -3.02 -3.02 -2.37 0.21 -2.74 -2.73 -1.98 -1.95 -2.04 -2.31 -2.31 -2.38
All values are relative to the QM/MM optimized crystal structure and are in kcal mol-1. For the crystal structure ligand pose, the electrostatic elec,)r elec,)4r interaction energies (eq 6), Eelec,)r ) -25.17, EQM/MM ) -29.02, Eelec,)r ) -3.85, Eelec,)4r ) -6.29, EQM/MM ) -7.25, and Eelec,)4r ) -0.96. gas pol gas pol a
parallels the QM/MM total energies. That is, the conformers with the lowest QM/MM energies are also among the lowest energy gas phase conformers. Somewhat unexpectedly, conformers 3, 4, and 12 and the crystal structure ligand conformer had the lowest gas-phase energies. The relative gas-phase energies are shown in Table 1. Conformers 3 and 4 were slightly more stable than the crystal structure (-0.87 and -0.88 kcal mol-1 lower in energy, respectively), and conformer 12 was higher in energy than the crystal structure conformer by 1.72 kcal mol-1. This result clearly shows that the preference for the binding mode observed in the crystal structure is largely determined by the conformational energetics (i.e., internal energy) of the ligand molecule. In other words, the other ligand conformers have simply adopted less favorable geometries. A comparison of the (unpolarized) gas-phase energies and the RMSD values for the 20 QM/MM optimized conformers is presented in Figure 3. The bottom panel shows the energetic trend in the gas phase conformers. It also includes the sum of the gas phase and van der Waals energies, and the middle panel shows the RMSD values for each conformer relative to the crystal conformer. Interestingly, this combination of terms largely captures the essential energetic aspects of binding. A gas-phase geometry optimization of the ligand was performed and the energy was compared to that of the ligand conformation in the optimized crystal structure. The crystal structure conformer was 12.3 kcal mol-1 higher in energy than the local minimum conformation obtained from the optimization. Although the conformation of the ligand in the crystal structure is significantly higher in energy than the gas-phase minimum energy conformer, it is relatively stable when compared to the other nonnative-like conformers. Trends in Energy Terms. In considering the individual energetic terms that comprise the QM/MM total energy, we see from Figure 3 that the polarized QM internal energy, E′, yields a trend similar to that of the gas-phase energy. Again, the three lowest energies correspond to conformers 3, 4, and 12. Additionally, a similar trend is observed for the van der Waals energies of the conformers. The lowest energy poses, conformers 3, 4, and 12, are sterically the most favorable ones. Curiously, there is an inverse relationship between the electrostatic interaction energies of the conformers and their corre-
sponding QM/MM total energies. Table 1 shows that conformers 3, 4, and 12 exhibit relatively small electrostatic stabilization energies, while those of the other conformers are much larger in magnitude and are therefore more favorable. The relatively unfavorable electrostatic interactions observed in the crystal conformer and the three lowest energy conformers do not significantly penalize them because their QM internal and van der Waals energies more than offset their relatively unfavorable electrostatics contributions in the QM/MM total energies. van der Waals Contributions to Ligand Binding. Significant stabilization of the crystallographic complex and conformers 3, 4, and 12 derives from the van der Waals interactions vdw for the between the protein and ligand. From Table 1, EQM/MM -1 crystal structure was -43.8 kcal mol , followed by conformers vdw 12, 3, and 4 with EQM/MM ) -38.8, -33.0, and -33.0 kcal -1 mol , respectively. Other conformers ranged from -17 to -25 kcal mol-1, indicating less favorable interactions with the protein. Role of Ligand Charge Polarization. Inclusion of ligand polarization has been cited as an important factor in accurately predicting docked structures.39,40 We set out to quantify the contribution of polarization associated with ligand binding for the current system under study. We computed the gas-phase energies and ESP charges of the 20 QM/MM-optimized ligand conformers and the optimized crystal structure at the B3LYP/ 6-31G* level. We then used the gas-phase ESP charges of each ligand conformer to compute the classical electrostatic interaction energies between the protein and ligand conformers. The difference between the energy computed with the ESP charges obtained in the protein environment and those obtained in the gas phase with the same geometries provides an approximation to the electrostatic polarization component of the QM/MM energy. The classical polarization energy is defined as Eelec pol ≡ elec(gas) elec(gas) elec,(r) - EQM/MM , where EQM/MM is the electrostatic interacE′QM/MM tion energy between the protein and ligand using gas-phase ESP charges for the ligand. This term was computed using distancedependent dielectric constants of ) rij and ) 4rij. The results for both values of are shown in Table 2. Eelec pol for the crystal structure conformer was -3.85 kcal mol-1 for ) rij. The only conformer with a less favorable
3174 J. Phys. Chem. B, Vol. 112, No. 10, 2008 value of Eelec pol was conformer 12, which had a polarization energy of -3.02, a difference of 0.83 kcal mol-1. The polarization energies for conformers 3 and 4 were lower than the crystal conformer by 0.56 kcal mol-1. All other ligand conformers were significantly more polarized with absolute values ranging from about -12 to -17.5 kcal mol-1. Our intention here was not to provide exact values for the polarization energies, but to show the overall trends with different values of the dielectric constant. These results suggest that, for the HCV polymerase-ligand system, the polarization energy is at a minimum, or a near-minimum, for the optimally bound pose. IV. Discussion Potential Energy Versus Free Energy. For the simulation of protein-ligand binding interactions, two primary options exist. One is to score and/or optimize ligand conformations on a potential energy surface. In this case, the ligand is flexible, and the protein may or may not be allowed to move. This approximation is the most common approach in the literature and is often sufficient for predicting structures of bound complexes. The other option is to allow the entire system to sample a free energy surface, typically with either explicit or continuum solvation. This option is accomplished by performing converged statistical mechanical sampling via molecular dynamics or Monte Carlo as in FEP calculations. This latter approach requires the use of a molecular mechanical, or possibly semiempirical QM representation of the ligand. Regardless, these methods are quite expensive computationally and obtaining converged results is difficult in practice. Additional errors may arise as a result of an inaccurate representation of the ligand. Thus, we have chosen to focus on potential energy surface optimization using an accurate quantum mechanical description of the ligand and a fixed protein structure. Our approximate treatment of solvation through the use of a distance-dependent dielectric constant provides satisfactory results in the present case, as evidenced by the agreement between our predicted binding pose and the crystal structure. Additionally, we would not be able to compare the energetics of various conformations directly if we allowed the protein to move. Ranking Conformers by QM/MM Total Energy. Other groups have used semiempirical QM/MM optimizations in conjunction with PB/SA calculations,41 or charges derived from high-level DFT calculations39 as input to provide improved electrostatics in their scoring functions. We have performed DFT-based QM/MM optimizations on a set of protein-ligand complexes and have used the QM/MM energies to rank the ligand conformers. It is possible that semiempirical methods could provide similar results with a reduced computational cost, although they are more likely to yield inaccurate energetics and structures. Simple Dielectric Models. Several of the conformers predicted by FlexX including conformer 1 only partially occupied the binding pocket. They instead adopted configurations that protruded into an open area that would otherwise be occupied by solvent. As a result, the carboxylate group of the ligand in these conformers preferred to interact strongly with Arg394, which is located outside of the binding pocket (see Figure 2). These predicted poses were problematic because they failed to account for strong protein-water and ligand-water interactions. Our QM/MM calculations would have suffered the same fate if we neglected the charge-charge damping effects provided by the distance-dependent dielectric constant. In this way, we were able to minimize the occurrence of false positives originating from long-range electrostatic interactions between
Parks et al. the ligand and charged hydrophilic residues at the proteinsolvent interface. The neglect or incomplete treatment of this effect may explain why FlexX predicted conformer 1 as the highest scoring binding pose. Another possibility for inaccurate scoring may be the tendency of the scoring function to favor hydrogen bonds between charged groups on the ligand and hydrogen bond donors on the protein. Protein-Ligand Electrostatic Interactions. Beierlein et al.42 found in their semiempirical QM/MM optimizations of a protein-ligand complex involving HIV-1 protease that the classical Coulomb interaction energy was -9.34 kcal mol-1 when polarized QM charges were used to compute the interactions. On the other hand, a value of -6.39 kcal mol-1 was obtained when the gas phase (i.e., unpolarized) QM charges were used. This implies an electrostatic polarization energy of -2.95 kcal mol-1. Our calculations yielded similar results, with an electrostatic polarization energy for the crystal structure of -3.85 kcal mol-1. While this result does demonstrate the effect of polarization, it does not show the interconnected nature of polarization and the quantum mechanical energy. For example, we see in Table 2 that many other conformers have significantly stronger polarization energies, but they also have higher QM/ MM total energies. This is because more favorable polarization energy comes at the expense of a higher internal energy and/or van der Waals energy. Thus, maximizing the polarization energy is not necessarily beneficial because there is a subtle interplay between other energetic terms. In other words, the deformation of the electron density arising from a charged environment results in a higher internal QM energy. Oxyanion Hole. HCV polymerase binds the ligand such that the carboxylate and one face of the adjacent benzamide group are solvent-exposed while the rest of the ligand occupies a hydrophobic region in the binding pocket (Figure 2). In the crystal structure, a water molecule forms a bridging hydrogen bond between one oxygen of the ligand carboxylate and the backbone N-H of Ser556. However, there are no direct hydrogen bonds between the carboxylate and the protein. In fact, the only direct hydrogen bond between the protein and the ligand involves the backbone N-H of Tyr448 and the carbonyl of the ligand benzamide group. This interaction stabilizes the orientation of the benzamide by helping to anchor it in place at the hydrophilic/hydrophobic interface. The hydrophilic interactions of the carboxylate and benzamide N-H group with the solvent remain favorable because there are no suitable hydrogen bond donors in the binding pocket. As a result, the carboxylate remains solvated, which precludes any energetic desolvation penalties that might otherwise occur. Based on the tendency of the carboxylate group to occupy this specific region, we term this site the oxyanion binding hole. Thus, the location of the oxyanion binding hole is influenced by the lack of suitable hydrogen bond donors within the binding site, the ability of the carboxylate and benzamide N-H groups to remain solvated and refrain from paying a desolvation penalty, and optimal electrostatic interactions of the charged ligand with the protein and solvent. Because this ligand has both hydrophilic and hydrophobic groups, related ligands that are more hydrophobic and lack polar fragments such as a carboxylate may not bind as effiently because of the proximity of the binding pocket to the solvent. V. Conclusions The crystallographic binding mode of the ligand (2Z)-2(benzoylamino)-3-[4-(2-bromophenoxy)phenyl]-2-propenoate in the binding pocket of HCV NS5B polymerase was analyzed
Hepatitis C Virus NS5B Polymerase via DFT-based QM/MM optimization and energetic decomposition. Our study reveals two key features. First, we successfully identified, from among the 20 initial docked poses, three poses that have low RMSD with respect to the crystallographic structure. Second, we were also able to decompose the energies of a number of protein-ligand complexes into contributions from the polarized QM internal energy of the ligand, QM/MM van der Waals, and QM/MM electrostatic energies in order to observe how the total energy and its components vary in optimal and sub-optimal geometric configurations of the ligands. Surprisingly, the internal energy determines the energy ordering of the docked poses, favoring three of the 20 docked poses. We then showed that solvation effects must be included to damp the electrostatic interactions properly and to obtain reliable, semiquantitative results. We also provide an explanation for the positioning of the oxyanion binding hole that is occupied by the carboxylate group of the ligand. This study demonstrates the validity of the QM/MM procedure in the HCV NS5B polymerase system and provides some understanding of the molecular basis for the observed docking specificity. Acknowledgment. The authors thank the National Institutes of Health (NIH R01-GM-061870) for financial support and Roche LLC for unrestricted support of the theory program in Chemistry at Duke University. References and Notes (1) Senn, H. M.; Thiel, W. QM/MM methods for biological systems. Top. Curr. Chem. 2007, 268, 173-290. (2) Warshel, A.; Levitt, M. Theoretical studies of enzymic reactions: Dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. J. Mol. Biol. 1976, 103, 227-249. (3) Zhang, Y.; Lee, T.-S.; Yang, W. A pseudobond approach to combining quantum mechanical and molecular mechanical methods. J. Chem. Phys. 1999, 110, 46-54. (4) Zhang, Y.; Liu, H.; Yang, W. Free energy calculation on enzyme reactions with an efficient iterative procedure to determine minimum energy paths on a combined ab initio QM/MM potential energy surface. J. Chem. Phys. 2000, 112, 3484-3492. (5) Cisneros, G. A.; Liu, H.; Zhang, Y.; Yang, W. Ab initio QM/MM study shows there is no general acid in the reaction catalyzed by 4-oxalocrotonate tautomerase. J. Am. Chem. Soc. 2003, 125, 10384-10393. (6) Cisneros, G. A.; Wang, M.; Silinski, P.; Fitzgerald, M. C.; Yang, W. The protein backbone makes important contributions to 4-oxalocrotonate tautomerase enzyme catalysis: Understanding from theory and experiment. Biochemistry 2004, 43, 6885-6892. (7) Liu, H.; Lu, Z.; Cisneros, G. A.; Yang, W. Parallel iterative reaction path optimization in quantum mechanical/molecular mechanical modeling of enzyme reactions. J. Chem. Phys. 2004, 121, 697-706. (8) Lu, Z.; Yang, W. Reaction path potential for complex systems derived from combined ab initio quantum mechanical and molecular mechanical calculations. J. Chem. Phys. 2004, 121, 89-100. (9) Hu, H.; Yang, W. Dual-topology/dual-coordinate free-energy simulation using a QM/MM force field. J. Chem. Phys. 2005, 123, 041102. (10) Hu, H.; Lu, Z.; Yang, W. QM/MM Minimum Free Energy Path: Methodology and application to triosephosphate isomerase. J. Chem. Theory Comput. 2007, 3, 390-406. (11) Becke, A. D. Density-functional thermochemistry. iii. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648-5652. (12) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. ReV. B 1988, 37, 785-789. (13) Grimme, S. Seemingly simple stereoelectronic effects in alkane isomers and the implications for Kohn-Sham density functional theory. Angew. Chem. Int. Ed. 2006, 45, 4460-4464. (14) Schreiner, P. R.; Fokin, A. A.; Pascal, R. A., Jr.; de Meijere, A. Many density functional theory approaches fail to give reliable large hydrocarbon isomer energy differences. Org. Lett. 2006, 8, 3635-3638. (15) Wodrich, M. D.; Corminboeuf, C.; Schreiner, P. R.; Fokin, A. A.; von Rague´ Schleyer, P. How accurate are DFT treaments of organic energies? Org. Lett. 2007, 9, 1851-1854. (16) Kristya´n, S.; Pulay, P. Can (semi)local density functional theory account for the London dispersion forces? Chem. Phys. Lett. 1994, 229, 175-180.
J. Phys. Chem. B, Vol. 112, No. 10, 2008 3175 (17) Zhao, Y.; Truhlar, D. G. Density functionals for noncovalent interaction energies of biological importance. J. Chem. Theory Comput. 2007, 3, 289-300. (18) Wu, Q.; Yang, W. Empirical correction to density functional theory for van der Waals interactions. J. Chem. Phys. 2002, 116, 515-524. (19) Jurecˇka, P.; C ˇ erny´, J.; Hobza, P.; Salahub, D. R. Density functional theory augmented with an empirical dispersion term. Interaction energies and geometries of 80 noncovalent complexes compared with ab initio quantum mechanics calculations. J. Comput. Chem. 2007, 28, 555-569. (20) Cui, Q.; Karplus, M. Quantum mechanical/molecular mechanical studies of the triosephosphate isomerase-catalyzed reaction: Verification of methodology and analysis of reaction mechanisms. J. Phys. Chem. B 2002, 106, 1768-1798. (21) Riccardi, D.; Li, G.; Cui, Q. Importance of van der Waals interactions in QM/MM simulations. J. Phys. Chem. B 2004, 108, 64676478. (22) Pfefferkorn, J. A.; Greene, M. L.; Nugent, R. A.; Gross, R. J.; Mitchell, M. A.; Finzel, B. C.; Harris, M. S.; Wells, P. A.; Shelly, J. A.; Anstadt, R. A.; Kilkuskie, R. E.; Kopta, L. A.; Schwende, F. J. Inhibitors of HCV NS5B polymerase. Part 1: Evaluation of the southern region of (2Z)-benzoylamino)-3-(5-phenyl-2-furyl)acrylic acid. Bioorg. Med. Chem. Lett. 2005, 15, 2481-2486. (23) Bohm, H. J. The development of a simple empirical scoring function to estimate the binding constant for a protein-ligand complex of known three-dimensional structure. J. Biol. Chem. 1994, 8, 243-256. (24) Kramer, B.; Rarey, M.; Lengauer, T. Evaluation of the FlexX incremental construction algorithm for protein-ligand docking. Proteins: Struct., Funct., Genet. 1999, 37, 228-241. (25) Rarey, M.; Kramer, B.; Lengauer, T.; Klebe, G. A fast flexible docking method using an incremented construction algorithm. J. Mol. Biol. 1996, 261, 470-489. (26) Vieth, M.; Hirst, J. D.; Kolinski, A.; Brooks, C. L., III Assessing energy functions for flexible docking. J. Comput. Chem. 1998, 19, 16121622. (27) Jain, T.; Cerutti, D. S.; McCammon, J. A. Configurational-bias sampling technique for predicting side-chain conformations in proteins. Protein Sci. 1998, 15, 2029-2039. (28) Huo, S.; Wang, J.; Cieplak, P.; Kollman, P. A.; Kuntz, I. D. Molecular dynamics and free energy analyses of cathepsin D-inhibitor interactions: Insight into structure-based ligand design. J. Med. Chem. 2002, 45, 1412-1419. (29) Ferrara, P.; Gohlke, H.; Price, D. J.; Klebe, G.; Brooks, C. L., III Assessing scoring functions for protein-ligand interactions. J. Med. Chem. 2004, 47, 3032-3047. (30) Gordon, J. C.; Myers, J. B.; Folta, T.; Shoja, V.; Heath, L. S.; Onufriev, A. H++: A server for estimating pKas and adding missing hydrogens to macromolecules Nucl. Acids Res. 2005, 33, W368-371. (31) Bashford, D.; Karplus, M. pKas of ionizable groups in proteins: Atomic detail from a continuum electrostatic model. Biochemistry 1990, 29, 10219-10225. (32) Myers, J.; Grothaus, G.; Narayana, S.; Onufriev, A. A simple clustering algorithm can be accurate enough for use in calculations of pKas in macromolecules. Proteins: Struct., Funct., Bioinf. 2006, 63, 928-938. (33) MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L., Jr.; Evanseek, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B., III, W. E. R.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wio´rkiewicz-Kuczera, J.; Yin, D.; Karplus, M. All-atom empirical potential for molecular modeling and dynamics studies of proteins. J. Phys. Chem. B 1998, 102, 3586-3616. (34) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general Amber force field. J. Comput. Chem. 2004, 25, 1157-1174. (35) Mann, G.; Yun, R. H.; Nyland, L.; Prins, J.; Board, J.; Hermans, J. The Sigma MD program and a generic interface applicable to multifunctional programs with complex, hierarchical command structure. In Computational Methods for Macromolecules: Challenges and Applications, Proceedings of the 3rd International Workshop on Algorithms for Macromolecular Modelling; Schlick, T., Gan, H. H., Eds.; Springer: New York, pp 129-145. (36) Hu, H.; Elstner, M.; Hermans, J. Comparison of a QM/MM force field and molecular mechanics force fields in simulations of alanine and glycine “dipeptides” (Ace-Ala-Nme and Ace-Gly-Nme) in water in relation to the problem of modeling the unfolded peptide backbone in solution. Proteins: Struct., Funct., Genet. 2003, 3, 451-463. (37) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; A.; M. J.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.;
3176 J. Phys. Chem. B, Vol. 112, No. 10, 2008 Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (38) Singh, U. C.; Kollman, P. A. An approach to computing electrostatic charges for molecules. J. Comput. Chem. 1984, 5, 129-145. (39) Cho, A. E.; Guallar, V.; Berne, B. J.; Friesner, R. A. Importance of accurate charges in molecular docking: Quantum mechanical/molecular mechanical approach. J. Comput. Chem. 2005, 26, 915-927.
Parks et al. (40) Friesner, R. A. Modeling polarization in proteins and protein-ligand complexes: Methods and preliminary results. AdV. Protein Chem. 2006, 72, 79-104. (41) Gra¨ter, F.; Schwarzl, S. M.; Dejaegere, A.; Fischer, S.; Smith, J. C. Protein/ligand binding free energies calculated with Quantum Mechanics/Molecular Mechanics. J. Phys. Chem. B 2005, 109, 1047410483. (42) Beierlein, F.; Lanig, H.; Schu¨rer, G.; Horn, A. H. C.; Clark, T. Quantum mechanical/molecular mechanical (QM/MM) docking: An evaluation for known test systems. Mol. Phys. 2003, 101, 24692480. (43) Humphrey, W.; Dalke, A.; Schulten, K. VMD - Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33-38. (44) Stone, J. An Efficient Library for Parallel Ray Tracing and Animation. Master’s thesis; University of Missouri: Rolla, 1998.
