MM Study

Binding of Benzylpenicillin to Metallo-β-lactamase: A QM/MM Study. Lars Olsen* ... Applications of First Principles Quantum Chemical Methods in Drug ...
11 downloads 0 Views 356KB Size
Download PDF
J. Phys. Chem. B 2004, 108, 17639-17648

17639

Binding of Benzylpenicillin to Metallo-β-lactamase: A QM/MM Study Lars Olsen,*,† T. Rasmussen,‡ L. Hemmingsen,§ and U. Ryde‡ Department of Medicinal Chemistry, The Danish UniVersity of Pharmaceutical Sciences, UniVersitetsparken 2, DK-2100 Copenhagen, Denmark, Department of Theoretical Chemistry, UniVersity of Lund, Chemical Center, P.O.B. 124 S-221 00 Lund, Sweden, and Department of Physics, The Quantum Protein Centre, The Technical UniVersity of Denmark, 2800 Kgs. Lyngby, Denmark ReceiVed: April 23, 2004; In Final Form: August 19, 2004

Metallo-β-lactamases are bacterial enzymes that may function with either one or two zinc ions bound in the active site. In this work, the binding of benzylpenicillin to mono-zinc metallo-β-lactamase from Bacillus cereus has been investigated in a docking procedure applying a combined quantum mechanical/molecular mechanical method as the final step. It is demonstrated that the substrate can bind with the carbonyl oxygen of the lactam ring coordinating to the zinc ion, and with the zinc-bound hydroxide ion in position for a nucleophilic attack on the carbonyl carbon of the lactam ring. In some structures, both the histidine and the cysteine at the other (unoccupied) metal-binding site are in a proper position to function as proton shuttles in proton transfer from the previously zinc-bound hydroxide, to the nitrogen in the lactam ring. In addition, the hydrophobic region formed by Phe34, Val39, Trp59, and Ala89 interacts with the phenyl group of benzylpenicillin, whereas the carboxylate group may be stabilized by Lys171 and Asn180. Alternatively, the carboxylate can bind to the zinc ion, prohibiting the nucleophilic attack of the zinc-bound hydroxide on the lactam carbonyl carbon. However, such a structure is energetically disfavored compared to the other enzymesubstrate complexes.

1. Introduction β-Lactamases are bacterial enzymes that catalyze the hydrolytic cleavage of the lactam C-N bond of β-lactam antibiotics. This inactivates the antibiotics, and the β-lactamases are one of the major reasons for bacterial resistance to these drugs. The β-lactamases have been divided into four classes, A-D.1 Enzymes belonging to class B depend on metal ions for their activity and pose a serious threat to public health, because they have a broad substrate profile and no inhibitors are clinically available.2 Metallo-β-lactamases have two potential Zn2+-binding sites. The enzymes from B. cereus (subclass B1) and A. hydrophila (subclass B2) bind the first Zn2+ ions with significantly higher affinity than the second Zn2+,3-7 whereas the enzymes from B. fragilis and S. maltophila bind them with similar binding affinity.8,9 The enzymes from B. cereus and B. fragilis are almost equally active with one or two Zn2+ ions,7,10,11 whereas the enzyme from S. maltophila is significantly more active and that from A. hydrophila is less active with the second Zn2+ bound.5,7,12 A number of crystal structures of metallo-β-lactamases have been published showing two metal ion binding sites, one of which must be occupied for activity.13-18 However, a major problem in the design of the inhibitors of metallo-β-lactamase is the lack of crystal structures with a native substrate bound. Therefore, interactions between metallo-β-lactamases and substrates or inhibitors have been studied with a variety of experimental and theoretical methods,13-31 see Table 1. In almost all of these studies, the β-lactam carbonyl oxygen, (OLac, see Figure 1) is found to interact with Zn1 and in most * Corresponding author. E-mail: [email protected]. † The Danish University of Pharmaceutical Sciences. ‡ University of Lund. § The Technical University of Denmark.

TABLE 1: Structural Studies of β-Lactamase: Substrate or Inhibitor Complexesa enzyme BCIIb BCII

substrate/inhibitors cephalosporin benzylpenicillin, meropenem benzylpenicillin cefotaxime thiomandelic acidc benzylpenicillin captopril

BCII BCII BCII BCIIb BCII, CphA CcrA ampicillin, ceftazidime, imipenem CcrA ampicillin, cephaloridine CcrAf biphenyl tetrazolec CcrA CcrA CcrA CcrA IMP-1 IMP-1

ampicillin mercaptocarboxylatec mercaptocarboxylatec captopril mercaptocarboxylatec succinic acid,c imipenem

method

ref

manual docking manual docking

15 14

MMd MD NMR docking, MD, QMd EXAFS

67 27 25, 26 22 29

manual docking

13

manual docking database screening, docking, X-ray diffraction docking, MD, QMd NMR MDe MM pol+CTf X-ray diffraction database screening, docking, MM, X-ray diffraction MD

30 23

IMP-1, cephalothin, cefotaxime, IMP-6 cephaloridine, ceftazidime L1 ampicillin, ceftazidime, manual docking, imipenem MC, MMd

19 20, 21 32 24 18 31 28 17

a CcrA, MBL from Bacteroides fragilis; BCII, MBL from Bacillus cereus; IMP-1, MBL from Pseudomonas aeruginosa; L1, MBL from Stenotrophomonas maltophilia. b Only one Zn2+ bound to the enzyme, otherwise two Zn2+ ions are bound. c Inhibitors. d Docking with the carboxylate of the thiazolidine ring bound to the Lys residue. e Two mutations, Ala171Thr and Asn208Asp, compared to CcrA genbank accession number M63556. f The SIBFA molecular mechanical method, which includes polarization and charge-transfer effects. MC ) Monte Carlo search for different conformers, MM ) molecular mechanical geometry optimization, MD ) molecular dynamics, QM ) quantum mechanical studies of the reaction.

10.1021/jp0482215 CCC: $27.50 © 2004 American Chemical Society Published on Web 10/15/2004

17640 J. Phys. Chem. B, Vol. 108, No. 45, 2004

Figure 1. Benzylpenicillin.

cases also with the side chain of the adjacent Asn residue (Asn180 in Bacillus cereus; following Fabiane et al.,14 we denote the zinc ion with three histidine ligands Zn1, and that with a histidine, an aspartate, and a cysteine ligand Zn2). Most of the suggested substrate-binding modes have an interaction between the β-lactam carboxylate moiety (COOPen) and a conserved Lys residue (Lys171 in B. cereus), and between the aromatic β-lactam side chain and a hydrophobic pocket (Phe34, Val39, Trp59, and Ala89 in B. cereus). Other possible protein-substrate contacts have been suggested; for example, COOPen may coordinate directly to Zn2, interact with His210, or interact with the so-called apical water molecule bound to Zn2.15,16,19 Moreover, it has been proposed that the nitrogen of the lactam ring, NLac, binds to Zn2,17 or that the ligand interacts with a flexible loop, residues 34-38 in B. cereus, that moves upon substrate binding.19-21,32 Previously, molecular dynamical (MD) simulations of the binding of benzylpenicillin22 and cefotaxime27 to mono-zinc β-lactamase from B. cereus have been performed. For benzylpenicillin, it was shown that OLac interacts with the backbone nitrogen of the conserved residue Asn180, whereas the side chain of the same residue interacts with the carbonyl group of the β-lactam 6-acylamino side chain, OPen. It was also found that COOPen of benzylpenicillin interacts with the side chain of Lys171, as well as with the zinc-bound hydroxide via a water molecule. Furthermore, hydrophobic contacts were made between the phenyl group of the antibiotics, PhePen, and the hydrophobic pocket in the protein that consists of Val39, Trp59, and Phe34. The distance between the zinc-bound OH- and the carbonyl carbon (CLac) of the lactam ring, on which a nucleophilic attack is performed, was relatively large, on average 4.9 Å, but in 2% of the sampled structures, a so-called near-attack conformer was found, with a distance below 4.0 Å. As compared to the binding of benzylpenicillin, cefotaxime binds not too far away from the zinc ion (4 Å), and the backbone of Asn180 is involved in the stabilization of the negatively charged COOgroup. In addition, CLac is found closer to zinc-bound OH- (3.2 Å) than for benzylpenicillin. Mutational analysis of the conserved Asn and Lys residues, which are often assumed to be involved in substrate binding, have been performed for the enzyme from Bacteriodes fragilis.33 The result indicated that the Lys residue plays a significant role in the binding of penicillin G, but it is not equally important in binding of cephalosporin, probably owing to differences in relative orientation of the carbonyl and carboxylate groups in the two substrates. The results for the Asn residue suggested a role in binding and activation of the substrate during catalysis. It is noteworthy that in this compilation of enzyme-ligand interaction studies, all except three15,22,27 have been devoted to the di-zinc form of the enzyme. Therefore, we here present a study of the binding of a typical substrate, benzylpenicillin, to the mono-zinc enzyme from B. cereus, employing a recent highresolution crystal structure. Moreover, it is apparent from the compilation that substrate binding may involve direct interactions with the zinc ion. All earlier docking studies have

Olsen et al. employed empirical (molecular mechanics) force fields, but it is well-known that it is hard to accurately describe the coordination of ligands to metal ions with such methods, because the results will depend strongly on the parameters used. This is particularly important for a zinc ion with its flexible coordination sphere, which may be 4-, 5-, or 6-coordinated with no preferred geometry. Therefore, we have employed accurate density functional calculations for this interaction in the final step of our docking studies by means of combined quantum mechanical and molecular mechanical (QM/MM) geometry optimizations, which at the same time take into account all other interactions of the substrate with the surrounding protein with molecular mechanics methods. Thus, the present results should give a more accurate picture of the binding of a penicillin substrate to a mono-zinc metallo-β-lactamase. With different starting points for the QM/MM geometry optimization, we obtain different binding modes, each of which provide information on possible enzyme-substrate interaction, which may serve to understand the function of the enzyme. 2. Methods Protein. Only one crystal structure of the metallo-β-lactamase from B. cereus has the flexible loop in the closed position, which is believed to be the case for the active enzyme. This is a binuclear Cd-structure at 1.35 Å resolution (1MQO).34 In this structure, a citrate molecule binds the two cadmium ions, while Asp90 unexpectedly coordinates to Cd1, probably as an effect of the binding of the citrate molecule to the active site. By superimposing 1MQO on a mono-zinc structure of metallo-βlactamase (1BVT from B. cereus at 1.85 Å16), we modified the position of Ala89 and Asp90 in 1MQO to overlay 1BVT, to obtain a starting structure with a closed loop. In addition, the cadmium ion bound to the three histidine residues was changed to zinc, the Zn1-bound water was extracted from 1BVT, and the second cadmium ion removed. This modified structure was the starting point in the docking study. On the basis of earlier studies, the following protonation states of Zn-ligands were selected: neutral, Hid84 (nitrogen in δ position protonated), Hie86 (nitrogen in  position protonated), Hid149, Cys168; negative, Asp90; positive, Hip210 (nitrogens in δ and  positions protonated).22 Upon inspection of the surroundings, the other three His residues in the protein were assigned as positive Hip residues. Docking with AutoDock. Benzylpenicillin was docked into the protein using AutoDock3.0.35 Amber36 charges were used for all atoms in the protein, except for the zinc ion, hydroxide, and the zinc-bound histidines. For these atoms, the charges determined by Diaz et al.22 were used. For benzylpenicillin, atomic charges were calculated at the B3LYP/6-31G* level. Default van der Waals parameters were used, except for the zinc ion, for which the parameters from Hoops et al.37 were used. For all other parameters, default values were used. Force-Field Calculations with Amber. Version 6 of Amber36 was used for the force-field calculations. No bonds were defined between the zinc ion and its ligands. In all energy minimizations, a convergence criterion of 0.001 kcal/mol on the energy and a maximum of 10 000 cycles was used. The MD simulations were carried out with fixed bond lengths using the SHAKE algorithm,38 and a time step of 2 fs. A cutoff for nonbonded interactions of 10 Å was used, and the nonbonded list was updated every 25 time steps. The temperature was 300 K and controlled using Berendsen’s coupling algorithm.39 An equilibration period of around 150 ps resulted in a stable trajectory, as judged from the potential energy of the system.

Binding of Benzylpenicillin to Metallo-β-lactamase

Figure 2. System 1. The QM system in the QM/MM calculations.

One of the enzyme-substrate complexes, denoted complex ADock, was selected from the AutoDock calculations. It was solvated in a water cap of 30 Å around the origin of the protein. This resulted in a complex consisting of 3543 protein atoms and 2489 water molecules. The positions of all water molecules and all hydrogen atoms of the complex were geometry optimized. This complex served as starting point for two MD simulations and will be called complex AMM. The input structure to the first MD simulation was obtained from complex AMM by geometry optimizing the atoms in benzylpenicillin, Zn2+, and OH-, subject to three restraints: OH--CLac, OPen to the backbone NH of Asp90, and OLacZn2+ using a flat-bottom potential with minimum/maximum distances of 1.7/2.3 Å. Subsequently, all atoms within a radius of 15 Å from Zn2+ were optimized without any restraints. By such a procedure, we avoided extensive changes of the active site during the initial steps of the MD simulation. Subsequently, an MD simulation of 1.5 ns was carried out without any restraints. A structure was geometry optimized after 750 ps; it is called complex BMD. The input structure of the second MD simulation was obtained from complex AMM by geometry optimizing all atoms in a radius of 15 Å from Zn2+, subject to the following restraints: Zn2+OLac, OLac-ND2 of Asn180, and COOPen-NZ of Lys171 (with the same potential as above). Then, an MD simulation of 750 ps was carried out with these restraints. Subsequently, a 1.2 ns MD simulation was carried out without any restraints. A structure was selected and geometry optimized (again all atoms within a radius of 15 Å from Zn2+), yielding complex CMD. Benzylpenicillin. For benzylpenicillin, force constants and charges were determined from DFT calculations using the B3LYP method40-42 with the 6-31G* basis set.43-47 To determine the force constants, the structure was optimized, the Hessian matrix was calculated, and the force constants were determined using the procedure of Seminario.48 RESP charges49 were fitted from the electrostatic potentials obtained by the Merz-Kollman scheme.50 The B3LYP calculations were carried out with Gaussian98.51 COMQUM. QM/MM geometry optimizations were carried out using the program ComQum.52-54 ComQum is a combination of a QM and an MM program. In this work, Turbomole55 was used as the QM program and Amber as the MM program. With ComQum, the molecular system is divided into three parts, systems 1, 2, and 3. System 1 is treated at the QM level, with the rest of the protein included as an array of point charges from the Amber library. System 1 is shown in Figure 2 and consisted of the imidazole rings of His84, His86, and His149, CH3COO- of Asp90, CH3SH of Cys168, Zn2+, and OH-. In addition, the lactam ring, the thiazolidine ring, the carboxylate, and the NH-CO-CH3 side chain of benzylpenicillin were

J. Phys. Chem. B, Vol. 108, No. 45, 2004 17641 included. This gives a total of 67 atoms in system 1. System 2 is optimized at the MM level and consisted of all residues within a radius of 8 Å of any atom in system 1. System 3 consisted of the rest of the protein and the 30-Å water sphere. It is treated at the MM level but is kept fixed during the optimizations. Special action is required for bonds between system 1 and 2 (the is latter denoted a junction atom). In the QM calculations on system 1, the junction atom is replaced by a hydrogen atom, the position of which is linearly related to that of the heavy junction atom of system 2. The following atoms were junction atoms: Cβ of His84, His86, and His149; CR of Asp90 and Cys168; the aliphatic carbon in the side chain of benzylpenicillin; and the methyl carbons attached to the five-membered ring of benzylpenicillin. The ratio of the distances C-Cjunction/ C-Hjunction was taken from the equilibrium bond lengths in the Amber force field and QM optimized structures of the amino acid fragment. Quantum chemical calculations were performed with the Becke-Perdew-86 functional56,57 along with DZpdf basis set on zinc58 and 6-31G*43-47 on the rest of the atoms. First, the geometry of system 1 was optimized, fully relaxing the geometry of system 2 in each cycle of the optimization. The convergence criteria were a change of 0.26 kJ/mol in the energy and 10-2 au in maximum norm of Cartesian gradients. Subsequently, system 2 was fixed, and system 1 was optimized until the changes were below 2.6 J/mol and 10-3 au. Poisson-Boltzmann Calculations. The interaction energy of system 1 with the surrounding protein and solvent was estimated by solving the Poisson-Boltzmann equation numerically using the MEAD 2.2 software.59,60 Charges for system 1 were obtained by fitting to the quantum electrostatic potential at points selected according to the Merz-Kollman (MK) scheme,50 whereas Amber36 partial charges were used for the remainder of the protein. The electrostatic potential for system 1 was obtained from a density perturbed by protein and explicit solvent point charges. Standard MK radii were used for fitting charges to the quantum system for all atoms except Zn2+, where a radius of 2.0 Å was used.61 The charge fitting calculations were performed with Gaussian9851 using B3LYP40-42 and the following combination of basis sets: 6-31G* for C and H, 6-31+G* for N, O, and S, as well as DZpdf for Zn.43-47,58 Parse radii were used for all atoms62 in the MEAD calculations. The system was split into three dielectric regions:  ) 1 for system 1,  ) 4 for the remainder of the protein, and  ) 80 for water. That is, system 1 is “solvated” in a dielectric continuum with two dielectric regions, one region with embedded charges modeling the protein and one region modeling the solvent. The chosen dielectric values have previously been shown to provide reasonable agreement with experimental results.63 Furthermore, both the electronic polarization and the nuclear rearrangement of system 1 are treated explicitly in the QM/MM calculations, and the polarization of system 1 is thus already reflected by the charges used for this dielectric region. Hence, no further polarization of system 1 is needed in the MEAD calculations, giving a natural choice of 1 for the dielectric value in this region. The MEAD calculations were based on the QM/MM optimized geometries, and the QM/MM protocol of replacing junction atoms with capping hydrogens was used. The partial charge of the atoms in the protein region that were bound to the junction atoms was set to zero, whereas the charge of the remaining atoms in the interface intersecting residues was adjusted so that the total charge of each of these residue moieties adds up to an integer value. This QM-MM intersection chargescreening protocol is also used in the QM part of the QM/MM

17642 J. Phys. Chem. B, Vol. 108, No. 45, 2004

Olsen et al.

TABLE 2: Contacts between Benzylpenicillin and the Proteina distance between heavy atoms (Å) Amber penicillin OLac

CLac COOPen OPen PhePen ThiazPen

protein Zn2+ OHAsn180(ND2) Asp90(NH) OHZn2+ Lys171(NZ) Asn180(N) Asp90(N) OHPhe34 Trp59 Ala89 Val39

AMM

BMD 2.30

QM/MM CMD

AQM/MM

BQM/MM

2.19

CQM/MM 2.50

2.92 3.08 2.55 3.02 4.28/4.29b 2.80 3.16 2.26 3.06 4.85 3.24 3.47

2.81 3.34 2.11

3.09

2.87

2.82/2.98 3.35

3.94/4.06b 2.76/3.90 3.00 3.16 3.98 4.99 3.35 4.36

2.94 3.33 5.26 3.45

3.52 3.82 3.28 4.55

3.43 2.31

2.99 2.73/3.05 2.89

3.41 5.36 3.51

3.56 3.70 3.50 4.26

a

Complexes AMM, BMD, and CMD originate from a docked structure (complex Adock), which is either geometry optimized (complex AMM) or from MD simulations (complexes BMD and CMD) with the Amber force-field. These three structures served as input for the QM/MM geometry optimizations carried out with Turbomole/Amber. For the definition of OLac, CLac, COOPen, and OPen, see Figure 1. PhePen and ThiazPen are short for the phenyl and thiazolidine rings of penicillin. For the polar interactions, a distance of 3.5 Å is used as the upper limit, and 5.5 Å is used for the hydrophobic interactions. For the hydrophobic interactions, the distance given in the table is the shortest one. b Water mediated interaction.

calculations.64 The screening was used both for the fitting of charges in system 1 and for the MEAD calculations. All explicit waters were removed in the MEAD calculations. The calculations were performed using 351 grid points with a spacing of 0.2 Å, centering the grid on the protein cap center. The reported solvation energies are averages of seven calculations in which the grid was moved 0.1 Å in both positive and negative direction along each Cartesian axis. The maximum difference among the seven calculations was less than 1 kcal/mol for each complex. The energies, that include solvation, presented in Table 3, (denoted EQM+PB with PB as an abbreviation for PoissonBoltzmann) were obtained by combining the solvation energies from MEAD, described above, with the B3LYP energies of system 1. The B3LYP energy calculations were performed with Gaussian9851 using the combination of basis sets described earlier in this section. We thus obtain a total energy that includes the energy of system 1 as well as the interaction energy between system 1 and the surroundings. However, the total energy does not include the energy of the protein part of the complex. Hence, differences in total energies do not reflect potential compensations by the protein of unfavorable interactions between system 1 and the surroundings. 3. Results In the following, we will describe the resulting structures of the various calculations. Section 3.1 describes the docking performed with AutoDock, resulting in an enzyme-substrate complex, complex ADock. As described in the Methods section, this structure was solvated, and all hydrogen atoms were optimized, yielding complex AMM. Complex AMM served as starting point for the MD simulations (yielding complex BMD, and CMD). Complexes AMM, BMD, and CMD will be described in Section 3.2. Subsequently, complexes AMM, BMD, and CMD were optimized with QM/MM methods, resulting in complexes AQM/MM, BQM/MM, and CQM/MM. These structures will be described in Section 3.3. 3.1. Docking of Benzylpenicillin into Metallo-β-lactamase. AutoDock was used to dock benzylpenicillin into metallo-βlactamase from B. cereus. A number of different conformations were obtained. The highest scoring structure had a CLac-OHdistance of ∼5 Å, COOPen bound to zinc, the methyl groups, which are attached to the thiazolidine ring, pointed toward the

Figure 3. One of the docked structures of benzylpenicillin using AutoDock, complex ADock. The arrow indicates the nucleophilic attack on CLac.

hydrophobic groove (Phe34, Val39, Trp59, and Ala89), and PhePen pointed toward the solvent. On the basis of the result of previous docking and mutation studies (Table 1), this was considered as a suboptimal binding mode and was therefore rejected. Thus, we used the second highest scoring docked structure, which has a short CLac-OH distance (around 3 Å). This conformation is shown in Figure 3. The phenyl ring of benzylpenicillin interacts with the hydrophobic region of metallo-β-lactamase, in which Phe34, Val39, Trp59, and Ala89 are found. In addition, OPen interacts with the backbone of Asp90, OLac with the side chain of Asn180, and COOPen with the backbone of Asn180. Lys171 points toward COOPen, but the distance is too long to form a real hydrogen bond. Interestingly, both Cys168 and His210 are relatively close to NLac. 3.2. Complexes AMM, BMD, and CMD. Complex ADock was used as starting point for the subsequent calculations. It was solvated, and the positions of all hydrogens and water molecules were optimized, resulting in complex AMM. This served as starting point for the MD simulations, from which complex BMD, and CMD, were obtained, and will be described below. Complex AMM. Details about some distances between the enzyme and the substrate for complex AMM are given in Table 2 and Figure 4. This structure is, of course, very similar to

Binding of Benzylpenicillin to Metallo-β-lactamase

J. Phys. Chem. B, Vol. 108, No. 45, 2004 17643

Figure 4. Contacts between substrate and enzyme in complex AMM. Complex AMM was obtained from complex ADock by adding water molecules and geometry optimizing their positions, as well as all other hydrogen atoms. Contacts between substrate and enzyme are indicated by the dashed lines.

complex ADock. Note, however, that, in complex AMM, there is a water-mediated interaction between Lys171 and COOPen. Complex BMD. Initial MD simulations of complex AMM showed large rearrangements in the structure around Zn2+ early in the simulation, probably owing to differences in the AutoDock and Amber force fields. Therefore, complex AMM was first geometry optimized using Amber, subject to three constraints (see Methods section) to maintain the conformation of penicillin. Subsequently, a 1.5 ns MD simulation was performed, without any restraints. In Figure 5, some important distances between enzyme and benzylpenicillin are shown as a function of time in the MD simulation. Zn2+ is coordinated by the three histidines and the hydroxide ion during the simulation. In addition, COOPen and OLac coordinate to Zn2+ during the simulations, and the zinc-bound OH- is found close to CLac (about 3.3 Å). Thus, the zinc ion is six-coordinated during the simulation. It can also be seen from Figure 5 that OPen interacts with the backbone N of Asp90, Cys168 is found near the hydroxide, and the initial hydrogen bond between His210 and OH- is broken shortly after the beginning of the simulation. In addition, the distance between Asp90 and OH- increases after around 750 ps, because a water molecule comes in as a bridge. PhePen is found in the hydrophobic groove formed by Phe34, Trp59, and Ala89. At the end of the MD simulation, large deviations in the position of PhePen are observed. This is due to an intramolecular rearrangement, in which PhePen changes position and interacts with the thiazolidine ring of penicillin. In the period 200-750 ps, energies and structures do not change much. A structure was subjected to energy minimization after 750 ps, and the resulting structure is denoted complex BMD. It is shown in Table 2 and Figure 6, and it shows the same interactions as described for the MD simulation. Complex CMD. As described in the Methods section, we carried out another MD simulation starting from complex AMM, using another set of restraints during the first 750 ns. This was done to keep some protein-ligand contacts before the simulation without restraints was run. In the following, we describe only the results of the unconstrained (1.2 ns) simulation. During this simulation, OLac coordinates to Zn2+, whereas COOPen forms a hydrogen bond to Lys171. Also, the backbone NH of Asn180 interacts with COOPen, but its position is more flexible. As observed in the MD simulation of complex B, PhePen

Figure 5. Selected contacts between the enzyme and penicillin (top: Polar interactions, middle: Hydrophobic interactions), and between the enzyme and the zinc-bound OH- (bottom) during the first MD simulation.

Figure 6. Contacts between substrate and enzyme in complex BMD.

interacts with Phe34, Trp59, Ala89, and less strongly with Val39. The structure around Zn2+ and OH- is similar to the one observed in the MD simulation of complex B, except that an additional water molecule coordinates to Zn2+ (distances not shown in Figure 7). Another notable difference is that His210 interacts with the zinc-bound OH-. Complex CMD was obtained by minimizing the structure obtained after 750 ps MD simulation, again because the energy and structure seemed to be stable at this point. Selected structural data and substrate-enzyme interactions are presented in Table

17644 J. Phys. Chem. B, Vol. 108, No. 45, 2004

Olsen et al.

Figure 9. Complex AQM/MM geometry optimized with QM/MM methods.

Figure 7. Selected contacts between the enzyme and penicillin (top: Polar interactions, middle: Hydrophobic interactions), and between the enzyme and the zinc-bound OH- (bottom) during the second MD simulation.

Figure 10. Complex BQM/MM geometry optimized with QM/MM methods.

it is in a perfect position for a nucleophilic attack on CLac (OHCLac ) 2.87 Å). Figure 8. Contacts between substrate and enzyme in complex CMD.

2 and Figure 8. The contacts are the same as those described during the simulation. 3.3. QM/MM Structures. Interactions between ligands coordinating to metal ions, such as Zn2+, are difficult to describe reliably using MM methods. Therefore, we used QM/MM methods (where this central part of the system is described by DFT) to optimize the structures obtained from MM, i.e., complexes AMM, BMD, and CMD. Complex AQM/MM. Complex AMM was geometry optimized at the QM/MM level, yielding complex AQM/MM. Selected structural data and substrate-enzyme interactions are presented in Figure 9 and Table 2. The overall binding mode of benzylpenicillin is similar to the structure obtained from the MM optimization. The most prominent difference is that the zinc-bound OH- now forms a hydrogen bond to OPen. Therefore,

Complex BQM/MM. Likewise, complex BMD was optimized at the QM/MM level, giving complex BQM/MM. Selected intermolecular distances of this complex are shown in Table 2 and Figure 10. The binding mode differs somewhat from the one obtained from the MD simulation: OLac does not coordinate to Zn2+ but instead points toward N in Asp90 (OLac-Asp90(N) ) 3.59 Å). This has the consequence that the hydrogen bond between this N atom and OPen is disrupted as compared to the MD structure (4.27 Å compared to 2.94 Å). OH- forms a hydrogen bond to OLac, which is probably the reason for the somewhat longer CLac-OH- (3.43 Å) compared to that of complex AQM/MM. To investigate the effect of another and more likely hydrogen bonding acceptor than OLac on the CLac-OHdistance, we also optimized another structure, in which the position of the hydrogen of OH- was modified to point in the direction of the carboxylate oxygen of Asp90, and reoptimized

Binding of Benzylpenicillin to Metallo-β-lactamase

J. Phys. Chem. B, Vol. 108, No. 45, 2004 17645 TABLE 3: Relative Energies for the Structures Geometry Optimized at the QM/MM Levela relative energy (kcal/mol) complex

∆EQM

∆EQM+Pol

∆EQM+PB

AQM/MM BQM/MM CQM/MM

3.9 1.6 0.0

0.9 6.5 0.0

12.8 26.2 0.0

a ∆E QM is the energy of system 1 determined at BP-86/6-31G* level with the DZpdf basis set on zinc. ∆EQM+Pol is the energy of system 1 including the surrounding point charges (BP-86/6-31G* level with the DZpdf basis set on zinc), where the point charge self energy is not included. ∆EQM+PB is the energy determined using Poisson-Boltzmann approach; see the Methods section for further details.

Figure 11. Complex CQM/MM geometry optimized with QM/MM methods.

the structure at the QM/MM level. This destabilized the complex by 5 kcal/mol, and the CLac-OH- distance increased to 4.08 Å. Complex CQM/MM. In complex CMD, an additional water molecule coordinated to Zn2+, as shown in Figure 8. This is probably an artifact of the MM force field parameters for zinc, and it has not been observed in any crystal structure. Therefore, it was removed. OH- was located almost exactly between the oxygen atom of Asp90 and the sulfur atom of Cys168. Its hydrogen atom pointed in the direction of Cys168, whose hydrogen, on the other hand, pointed toward its own backbone carbonyl oxygen. Interestingly, Cys168 was relatively close to both COOPen and NLac and could in principle form hydrogen bonds to either of these two groups with its HG atom. Therefore, both starting points were tried, but neither of them converged. Since the OH- is positioned between the side chain of Asp90 and Cys168, we also tried to start the QM/MM optimization with OH- forming a hydrogen bond to Asp90. With this starting point, the optimization converged. Three different starting points were used, differing in whether the HG atom of Cys168 formed a hydrogen bond to NLac, COOPen, or OH-. The optimized structure for the first two starting points ended up with similar structures, with significantly weakened hydrogen bonds. For the third starting point, the hydrogen bond between HG and OHwas broken (probably because the hydrogen on NE2 of His210 already forms a hydrogen bond to OH-), and in the optimized structure, the position of the HG atom was is not too far away from NLac and COOPen, i.e., similar to the structures obtained from the two other starting points. In the following, the structure from the first starting point will be described, because it was lowest in energy (Figure 11). The structure is similar to the one obtained from the MD simulation. The Cys168(HG)-NLac distance of 2.70 Å is too long to be a real hydrogen bond, but it is in a good position to donate a proton to NLac. As in complex A, OH- is close to CLac (2.99 Å). Energies. We have now obtained several different possible structures of benzylpenicillin docked to metallo-β-lactamase. To obtain information about the relative stability of the various structures, we have studied their energies. This is far from a trivial task, because the use of MD simulations and subsequent QM/MM geometry optimizations cause large differences in the protein structure, and estimates should ideally involve entropic effects (free energies).

Energies from the QM/MM calculations are presented in Table 3. The first column shows the quantum chemical energy of the quantum system (∆EQM), evaluated in a vacuum. It gives a direct and accurate estimate of the relative stability of the quantum system (including the zinc coordination) in the various complexes. It indicates that complex CQM/MM is slightly more stable than the other two complexes (by 2-4 kcal/mol). These energies do not contain any information about the influence of the surrounding on the active site. The next column in Table 3 contains the quantum chemical energy of the quantum system, calculated in the field of the surrounding protein, modeled as an array of point charges (∆EQM+Pol). The self-energy of the point charges are not included. It can be seen that the surrounding protein stabilizes complex AQM/MM relative the other complexes, but that complex CQM/MM is still the most stable one (by 1-2 kcal/mol). These energies indicate how the surrounding protein stabilizes the various complexes by electrostatic interactions. However, they give no information about the stability of the surrounding protein itself for the three complexes. This energy is much harder to estimate, partly because of the differences in protein structure, but also because the water molecules inside the protein may adopt different positions. Of course, the ∆EQM+Pol is subject to the same problem, but there the effect is second-order and therefore has a relatively small influence. We used the Poisson-Boltzmann approach to consider the rest of the protein in the calculations. In this approach, the structure of the surrounding protein is taken from the QM/MM calculation, but solvent molecules are replaced by a continuum solvent. This way, we automatically obtain an ensemble average for the free energy of the solvent. In addition, the electrostatic interactions between the quantum system and the protein are damped by a dielectric constant of 4. Energies obtained in this way are also included in Table 3 (∆EQM+PB) and show that complex CQM/MM is strongly stabilized by the protein. The results indicate that it is 13 kcal/mol more stable than complex AQM/MM and 26 kcal/mol more stable than complex BQM/MM. These energies are our best estimates of the relative stabilities of the various complexes. Thus, all energy estimates point out the same complex as the most stable, complex CQM/MM, which gives some confidence to this conclusion. 4. Discussion The study of binding of substrates to metallo-β-lactamases at the atomic level is a central issue, both as a starting point for understanding and modeling of the reaction catalyzed by these enzymes, and to gain insight into potential enzyme-inhibitor interactions to be used in structure-based design of inhibitors. The compilation in the Introduction has indicated a number of important interactions between metallo-β-lactamase and the substrates (cf. Table 1). The aim of this investigation is to find

17646 J. Phys. Chem. B, Vol. 108, No. 45, 2004 structures of benzylpenicillin bound to mono-zinc metallo-βlactamase from B. cereus, which fulfill many of the suggested interactions and compare their respective energies. A particularly important structural feature of the complexes is the distance between CLac and the zinc-bound hydroxide ion, because the latter is assumed to perform the nucleophilic attack as the first step of the reaction leading to hydrolysis of the lactam ring. Thus, we wish to investigate if this distance is sufficiently short for this reaction path to be feasible. The docking has been performed using a sequence of three methods with increased sophistication: docking with a simple energy function, MD simulation, and QM/MM calculations. In the latter two methods, it is unlikely that major changes in the binding modes take place. Therefore, we have tested various binding modes and sets of interactions by introducing restraints early in the MD simulations or by explicitly changing the hydrogen-bond patterns before the QM/MM minimizations. The systems used as starting points for the QM/MM geometry optimizations were selected from the MD simulation after 750 ps because the structures were stable in this time regime, and thus, the selected structures are representative for a large fraction of the structures sampled during the MD run. This strategy was chosen because computational ressources dictate that only a few structures were selected for the QM/MM geometry optimizations. Ideally, an exhaustive sampling should be carried out, also at the QM/MM level. We have obtained three possible structures for the binding of benzylpenicillin to mono-zinc metallo-β-lactamase from B. cereus. All structures have in common the binding of PhePen in the hydrophobic groove of the protein (involving Phe34, Trp59, and Ala87), see Figures 9-11. However, the binding modes of the substrate differ in several ways. (1) The zinc ion: In complex BQM/MM, it is coordinated by COOPen, and in complex CQM/MM it is coordinated by OLac, whereas in complex AQM/MM, the substrate does not coordinate at all to the metal ion. (2) OLac: In complex AQM/MM, it forms hydrogen bonds to Asn180(HD2), and in complex BQM/MM it forms a hydrogen bond to OH-, whereas in complex CQM/MM, it binds to Zn2+. (3) COOPen: In complex AQM/MM and complex CQM/MM, it interacts with the side chain of Lys171 (a water-mediated interaction in the former complex) and the backbone of Asn180, whereas in complex BQM/MM, it binds to Zn2+. (4) NLac: In complex CQM/MM, it is close to Cys168, whereas in complex AQM/MM both Cys168 and His210 are sufficiently close to be considered as potential proton shuttles. Protonation of NLac is believed to be the next step, after the nucleophilic attack on CLac in the lactam hydrolysis. For complex AQM/MM and complex CQM/MM, the zinc-bound hydroxide seems perfectly positioned for a nucleophilic attack on CLac, and the reaction may proceed along the path on which there is currently consensus.22,65,66 Complex CQM/MM has the lowest energy (Table 3) and reproduces several enzymesubstrate contacts proposed in the literature: OLac coordinates to the zinc ion, and the zinc-bound hydroxide seems perfectly positioned for nucleophilic attack on CLac. The distance between the zinc-bound hydroxide and CLac is 2.99 Å, which is even closer that the gas-phase optimum distance of 3.21 Å.65 The two conserved residues Lys171 and Asn180 participate in the enzyme-substrate interaction, which emphasizes the importance of these two residues in the binding of the substrate. The notion that Asn180 is also involved in the activation of the substrate through interaction with OLac33 is supported by complex AQM/MM, but it is not found in complex CQM/MM. In addition, Cys168 and His210 are in positions that allow them to act as proton

Olsen et al. shuttles in the proton transfer to NLac, which occurs after the nucleophilic attack. Finally, the hydrophobic interactions described above, common for all three structures, are probably central in achieving a high affinity between substrate and enzyme. Complex BQM/MM, however, is very different from the two other complexes, and also from what has been previously observed in the literature. It does not seem optimal for the enzyme to accommodate the substrate in this binding mode, because the CLac-OH- distance is longer, and therefore, the reaction is less likely to occur. It is not straightforward to visualize reaction paths where this binding mode is productive, but one more complex mechanism could involve jumping of the metal ion between the zinc ion binding sites. Three other studies have been devoted to the binding of substrate to mono-zinc metallo-β-lactamase.15,22,27 Carfi et al.15 studied the binding of a cephalosporin by manually docking. Their results suggested that OLac binds to the zinc ion in agreement with complex CQM/MM. They also suggest that COOPen is stabilized by His210, an interaction not found in any of our optimizations of benzylpenicillin. However, this is probably because COOPen is positioned differently in cephalosporins. It is somewhat surprising that our structures obtained during the MD simulations differ from those presented by Diaz et al.,22 as the same force field, Amber, was used in both studies. In their study, the substrate did not bind as close to the zinc ion as found in this work, and only 2% of their structures has the CLac-OH- distance of less than 4 Å. There may be two reasons for this: One is that different enzymes were studied, and the other is that the parametrization of benzylpenicillin differs. Considering that the QM/MM method should be more accurate than MM and that our structures seem more reasonable from a mechanistic point of view, we believe that the structures presented in this work give a better description of the binding of the substrate to mono-zinc metallo-β-lactamase. Finally, as compared to the binding of cefotaxime to B. cereus carried by Dal Peraro et al.,27 complex CQM/MM has several similarities, e.g., relatively short CLac-OH- as well as Lys171 and Asn180 interactions with COOPen. However, Zn2+ binds to OLac with a somewhat longer distance in their study (around 4 Å) than it does for benzylpenicillin. We believe that the binding modes in complex CQM/MM and complex AQM/MM are more likely candidates for the correct binding mode of benzylpenicillin to metallo-β-lactamase than complex BQM/MM. This is confirmed by the energies in Table 3, which show that complex BQM/MM is disfavored relative to complex AQM/MM and complex CQM/MM. According to these energies, the latter complex should be the best candidate for the binding mode of the substrate in metallo-β-lactamase. However, as the quantum chemical energies are quite similar, see Table 3, it cannot be excluded that there are several binding modes of the substrate, as has been repeatedly been suggested in the literature.17,24,32 5. Conclusions We have performed a docking of benzylpenicillin to the active site of mono-zinc metallo-β-lactamase from B. cereus. The results indicate that the best binding mode has OLac coordinated to the zinc ion, that COOPen forms hydrogen bonds to Lys171 and Asn180, and that PhePen binds in the hydrophobic groove. This is in accordance with the fact that Lys171 and Asn180 are conserved in all known metallo-β-lactamases and with mutation studies, showing important effects of these two residues. It is also in accordance with most earlier investigations of the protein (cf. Table 1). In this binding mode, CLac is positioned close

Binding of Benzylpenicillin to Metallo-β-lactamase (2.99 Å) to the zinc-bound hydroxide ion, which is believed to perform the nucleophilic attack on the substrate. This is even closer than the gas-phase optimum distance of the reactant complex (3.21 Å).65 Moreover, Cys168 and His210 are in ideal positions to work as proton shuttles during the breakdown of the tetrahedral intermediate. This lends strong credence to this docking mode. We currently try to follow the reaction mechanism in the protein using similar QM/MM methods. We would also like to point out the importance of using QM/ MM methods in the docking. We have seen that the MM and MD structures differ quite appreciably from the QM/MM structures, especially for the interaction with the zinc ion. In particular, there seems to be a strong tendency to a high coordination number (up to six) of Zn2+ in the classical simulations. Moreover, we have seen that QM/MM energies are quite different from the MM energies, especially for the active site. Altogether, these results show the importance of using accurate methods to describe interactions around metal ions. Acknowledgment. These studies have been supported by computer resources of Lunarc at Lund University. L.O. acknowledges financial support from Novo Nordisk A/S (Novo Nordisk prize). T.R. acknowledges financial support from the Carlsberg Foundation. U.R. acknowledges financial support from the Swedish research council. References and Notes (1) Galleni, M.; Lamotte-Brasseur, J.; Rossolini, G. M.; Spencer, J.; Dideberg, O.; Fre`re, J.-M. Antimicrob. Agents Chemother. 2001, 45, 660. (2) Matagne, A.; Dubus, A.; Galleni, M.; Fre`re, J. M. Nat. Prod. Rep. 1999, 16, 1. (3) Orellano, E. G.; Girardini, J. E.; Cricco, J. A.; Ceccarelli, E. A.; Vila, A. J. Biochemistry 1998, 37, 10173. (4) Paul-Soto, R.; Bauer, R.; Fre´re, J.-M.; Galleni, M.; Meyer-Klaucke, W.; Nolting, H.; Rossolini, G. M.; de Seny, D.; Hernandez-Valladares, M.; Zeppezauer, M.; Adolph, H.-W. J. Biol. Chem. 1999, 274, 13242. (5) Valladares, M. H.; Kiefer, M.; Heinz, U.; Paul-Soto, R.; MeyerKlaucke, W.; Nolting, H. F.; Zeppezauer, M.; Galleni, M.; Fre`re, J.-M.; Rossolini, G. M.; Amicosante, G.; Adolph, H.-W. FEBS Lett. 2000, 467, 221. (6) de Seny, D.; Heinz, U.; Wommer, S.; Kiefer, M.; Meyer-Klaucke, W.; Fre`re, M.; Galleni, J.-M; Bauer, R.; Adolph, H.-W. J. Biol. Chem. 2001, 276, 45065. (7) Wommer, S.; Rival, S.; Heinz, U.; Galleni, M.; Fre`re, J.-M.; Franceschini, N.; Amicosante, G.; Rasmussen, B.; Bauer, R.; Adolph, H.W. J. Biol. Chem. 2002, 277, 24142. (8) Crowder, M. W.; Wang, Z.; Franklin, S. L.; Zovinka, E. P.; Benkovic, S. J. Biochemistry 1996, 35, 12126. (9) Crowder, M. W.; Walsh, T. R.; Banovic, L.; Petit, M.; Spencer, J. Antimicrob. Agents Chemother. 1998, 42, 921. (10) Rasia, R. M.; Vila, A. J. Biochemistry 2002, 41, 1853. (11) Paul-Soto, R.; Hernandez-Valladares, M.; Galleni, M.; Bauer, R.; Zeppezauer, M.; Fre`re, J.-M.; Adolph, H.-W. FEBS Lett. 1998, 438, 137. (12) Hernandez-Valladares, M.; Felici, A.; Weber, G.; Adolph, H.-W.; Zeppezauer, M.; Rossolini, G. M.; Amicosante, G.; Fre`re, J.-M.; Galleni, M. Biochemistry 1997, 36, 11534. (13) Concha, N. O.; Rasmussen, B. A.; Bush, K.; Herzberg, O. Structure 1996, 4, 823. (14) Fabiane, S. M.; Sohi, M. K.; Wan, T.; Payne, D. J.; Bateson, J. H.; Mitchell, T.; Sutton, B. J. Biochemistry 1998, 37, 12404. (15) Carfi, A.; Pares, S.; Due´e, E.; Galleni, M.; Duez, C.; Fre`re, J.-F.; Dideberg, O. EMBO J. 1995, 14, 4914. (16) Carfi, A.; Duee, E.; Galleni, M.; Fre`re, J.-M.; Dideberg, O. Acta Crystallogr. 1998, D54, 313. (17) Ullah, J. H.; Walsh, T. R.; Taylor, I. A.; Emery, D. C.; Verma, C. S.; Gamblin, S. J.; Spencer, J. J. Mol. Biol. 1998, 284, 125. (18) Concha, N. O.; Janson, C. A.; Rowling, P.; Pearson, S.; Cheever, C. A.; Clarke, B. P.; Lewis, C.; Galleni, M.; Fre`re, J.-M.; Payne, D. J.; Bateson, J. H.; Abdel-Meguid, S. S. Biochemistry 2001, 39, 4288. (19) Krauss, M.; Gresh, N.; Antony, J. J. Phys. Chem. B 2003, 107, 1215. (20) Scrofani, S. D. B.; Chung, J.; Huntley, J. J. A.; Benkovic, S. J.; Wright, P. E.; Dyson, H. J. Biochemistry 1999, 38, 14507.

J. Phys. Chem. B, Vol. 108, No. 45, 2004 17647 (21) Huntley, J. J. A.; Scrofani, S. D. B.; Osborne, M. J.; Wright, P. E.; Dyson, H. J. Biochemistry 2000, 39, 13356. (22) Dı´az, N.; Sua´rez, D.; Merz, K. M., Jr. J. Am. Chem. Soc. 2001, 123, 9867. (23) Toney, J. H.; Fitzgerald, P. M. D.; Olson, N.; Grover-Sharma, S. H.; May, W. J.; Sundelof, J. G.; Vanderwall, D. E. Chem. Biol. 1998, 5, 185. (24) Antony, J.; Gresh, N.; Olsen, L.; Hemmingsen, L.; Schofield, C. J.; Bauer, R. J. Comput. Chem. 2002, 23, 1281. (25) Mollard, C.; Moali, C.; Papamicael, C.; Damblon, C.; Vessilier, S.; Amicosante, G.; Schofield, C. J.; Galleni, M.; Fre`re, J.-M.; Roberts, G. C. K. J. Biol. Chem. 2001, 276, 45015. (26) Damblon, C.; Jensen, M.; Ababou, A.; Barsukov, I.; Papamicael, C.; Schofield, C. J.; Olsen, L.; Bauer, R.; Roberts, G. C. K. J. Biol. Chem. 2003, 278, 29240. (27) Dal Peraro, M.; Vila, A. J.; Carloni, P. Proteins: Struct., Funct., Bioinf. 2004, 54, 412. (28) Oelschlaeger, P.; Schmid, R. D.; Pleiss, J. Biochemistry 2003, 42, 8945. (29) Heinz, U.; Bauer, R.; Wommer, S.; Meyer-Klaucke, W.; Papamichaels, C.; Bateson, J.; Adolph, H. W. J. Biol. Chem. 2003, 278, 20659. (30) Carfi, A.; Duee, E.; Paul-Soto, R.; Galleni, M.; Fre`re, J.-M.; Dideberg, O. Acta Crystallogr. 1998, D54, 47. (31) Toney, J. H.; Hammond, G. G.; Fitzgerald, P. M. D.; Sharma, N.; Balkovec, J. M.; Rouen, G. P.; Olson, S. H.; Hammond, M. L.; Greenlee, M. L.; Gao, Y. D. J. Biol. Chem. 2001, 276, 31913. (32) Salsbury, F. R., Jr.; Crowley, M. F.; Brooks, C. L., III. Proteins: Struct., Funct., Genet. 2001, 44, 448. (33) Yanchak, M. P.; Taylor, R. A.; Crowder, M. W. Biochemistry 2000, 39, 11330. (34) Garcia-Saez, I.; Chantalat, L.; Dideberg, O. Manuscript in preperation. (35) Morris, G. M.; Goodsell, D. S.; Halliday, R. S.; Huey, R.; Hart, W. E.; Belew, R. K.; Olson, A. J. J. Comput. Chem. 1998, 19, 1639. (36) Case, D. A.; Pearlman, D. A.; Caldwell, J. W.; Cheatham, T. E., III; Ross, W. S.; Simmerling, C. L.; Darden, T. A.; Merz, K. M., Jr.; Stanton, R. V.; Cheng, A. L.; Vincent, J. J.; Crownley, M.; Tsui, V.; Radmer, R. J.; Duan, Y.; Petera, J.; Massova, I.; Seibel, G. L.; Singh, U. C.; Weiner, P. K.; Kollman, P. A. AMBER 6; University of California: San Francisco, 1999. (37) Hoops, S. C.; Anderson, K. W.; Merz, K. M., Jr. J. Am. Chem. Soc. 1991, 113, 8262. (38) Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (39) Berendsen, H. J. C.; Postma, H. P. M.; van Gunsteren, W. F.; Di Nola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (40) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. 1988, B37, 785. (41) Becke, A. D. J. Chem. Phys. 1993, 98, 1372. (42) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (43) Hehre, W. J.; Ditchfield, R.; Pople, J. A. J. Chem. Phys. 1972, 56, 2257. (44) Ditchfield, R.; Hehre, W. J.; Pople, J. A. J. Chem. Phys. 1971, 54, 724. (45) Hariharan, P. C.; Pople, J. A. Mol. Phys. 1974, 27, 209. (46) Gordon, M. S. Chem. Phys. Lett. 1980, 76, 163. (47) Hariharan, P. C.; Pople, J. A. Theor. Chim. Acta 1973, 28, 213. (48) Seminario, J. M. Int. J. Quantum Chem. 1996, 60, 59. (49) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269. (50) Besler, B. H., Jr.; Merz, K. M., Jr.; Kollman, P. A. J. Comput. Chem. 1990, 11, 431. (51) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A., Jr.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B. G.; Chen, W.; Wong, M. W.; Andres, J. L.; Head-Gordon, M.; Replogle, E. S.; Pople, J. A. Gaussian 98, revision A.5; Gaussian, Inc.: Pittsburgh, PA, 1998. (52) Ryde, U. J. Comput.-Aided Mol. Des. 1996, 10, 153. (53) Ryde, U.; Olsson, M. H. M. Int. J. Quantum Chem. 2001, 81, 335. (54) Ryde, U.; Olsen, L.; Nillson, K. J. Comput. Chem. 2002, 23, 1058. (55) Treutler, D.; Ahlrichs, R. J. Chem. Phys. 1995, 102, 346. (56) Perdew, J. P. Phys. ReV. B 1986, 33, 8822. (57) Becke, A. D. Phys. ReV. A 1988, 38, 3098. (58) Scha¨fer, A.; Horn, H.; Ahlrich, R. J. Chem. Phys. 1992, 97, 2571.

17648 J. Phys. Chem. B, Vol. 108, No. 45, 2004 (59) Bashford, D.; Gerwert, K. J. Mol. Biol. 1992, 224, 473. (60) Bashford, D. Scientific Computing in Object-Oriented Parallel EnVironments; Springer: Berlin, 1997; Vol. 1343. (61) Sigfridsson, E.; Ryde, U. J. Comput. Chem. 1998, 19, 377. (62) Sitkoff, D.; Sharp, K. A.; Honig, B. J. Phys. Chem. 1994, 98, 1978. (63) Bashford, D.; Case, D. A.; Dalvit, C.; Tennant, L.; Wright, P. E. Biochemistry 1993, 32, 8045.

Olsen et al. (64) Ryde, U.; Hemmingsen, L. J. Biol. Inorg. Chem. 1997, 2, 567. (65) Olsen, L.; Antony, J.; Ryde, U.; Adolph, H.-W.; Hemmingsen, L. J. Phys. Chem. B 2003, 107, 2366. (66) Wang, Z.; Fast, W.; Valentine, A. M.; Benkovic, S. J. J. Curr. Opin. Chem. Biol. 1999, 3, 614. (67) Prosperi-Meys, C.; Wouters, J.; Galleni, M.; Lamotte-Brasseur, J. Cell. Mol. Life Sci. 2001, 58, 2136.