Supplementary Material (ESI) for PCCP This journal is © The Owner Societies 2000 ================================================== Molecular dynamics study of aromatic oxygenation reactions catalysed by the enzyme p-hydroxybenzoate hydroxylase Salomon R. Billeter (a), Christof F. W. Hanser (a), Tiziana Z. Mordasini (a), Mirjam Scholten (b), Walter Thiel* (b) and Wilfred F. van Gunsteren (c) a Organisch-chemisches Institut, Universitaet Zuerich, Winterthurerstrasse 190, CH-8057 Zuerich, Switzerland b Max-Planck-Institut fuer Kohlenforschung, Kaiser-Wilhelm-Platz 1, D-45470 Muelheim an der Ruhr, Germany c Laboratorium fuer Physikalische Chemie, ETH Zuerich, ETH-Zentrum, CH-8092 Zuerich, Switzerland Contents -------- 1. Validation of AM1: Comparison with ab initio and DFT results 2. Unrestrained environment in MD simulations 3. Definition of the MM' interaction region in MD simulations 4. Technical details of QM/MM geometry optimisations in PHBH 5. Technical note concerning Figure 5 of the paper 1. Validation of AM1: Comparison with ab initio and DFT results ---------------------------------------------------------------- We have considered the substrates p-oxybenzoate dianion (1), p-hydroxybenzoate anion (2), and p-oxybenzoic acid anion (3), and their reactions with hydrogen peroxide as discussed in section IV. The formulas for substrates (1)-(3) and for the hydroxylation products (1a)-(3a) are shown in Figure 1 of the paper. The alternative oxygenation products (1b)-(3b) are derived from (1a)-(3a) by abstracting a proton from the OH group at C3. The following reactions have been studied: (1) + HOOH -> (1a) + OH- (analogously for (2) and (3)) (1) + HOOH -> (1b) + HOH (analogously for (2) and (3)) The validation studies employed ab initio Hartree-Fock (HF) [S1] and hybrid density functional (B3LYP) [S2] calculations with standard basis sets [S1]. The 6-31++G(d,p) basis is derived from the 6-31G(d) basis by adding two sets of diffuse functions for a proper description of anions and one set of p-polarisation functions at hydrogen. The results with the larger 6-31++G(d,p) basis are expected to be more reliable. All HF and B3LYP calculations were performed using the Gaussian98 program [S3] and involved full geometry optimisations. Table S1 lists the relative energies for the model system (1)+HOOH, employing the same notation as in section IV of the paper. Table S1. Relative energies (kcal/mol) for the reactants (R-sep, R-compl), the transition state (TS), and the products (P-compl, P-sep) in the reaction between the p-OHB dianion (1) and hydrogen peroxide. R-sep R-compl TS P-sep P-compl P-sep (1) (1) (1a) (1b) (1b) HF/6-31G(d) 0.0 -17.8 38.2 -62.8 -62.1 -32.8 HF/6-31++G(d,p) 0.0 -28.7 37.1 -78.8 -62.0 -37.0 B3LYP/6-31G(d) 0.0 -21.5 -1.3 -50.3 -59.0 -26.3 B3LYP/6-31++G(d,p) 0.0 -16.4 -1.7 -73.0 -58.5 -31.2 AM1 0.0 -18.1 21.6 -60.3 -55.8 -37.3 The reference results from the HF and B3LYP calculations show some scatter and do not yet provide results that are converged with regard to basis set extension and electron correlation. The AM1 results lie generally in the same range as the reference results and appear to be generally satisfactory. The largest uncertainty is associated with the energy of TS where the HF and B3LYP results differ by almost 40 kcal/mol, with AM1 approximately halfway in between. It is known that HF tends to overestimate barriers while B3LYP is generally more realistic in this regard. To check this further, we have computed single-point MP2/6-31++G(d,p)//HF/6-31++G(d,p) energies and find that TS lies 27.2 kcal/mol above R-compl at this level (compared with 65.8 kcal/mol for HF/6-31++G(d,p), 14.7 kcal/mol for B3LYP/6-31++G(d,p), and 39.7 kcal/mol for AM1). Hence, both ab initio HF and AM1 overestimate this barrier compared with the correlated MP2 and B3LYP methods. Our validation attempts for substrates (2) and (3) are far less extensive. Table S2 contains the corresponding relative energies. Table S2. Relative energies (kcal/mol) involving the substrates (2) and (3) using an analogous notation as in Table 1. Reference system (3) R-sep(2) R-sep(3) R-sep(3) R-sep(3) Relative energy for (2) P-sep(2a) P-sep(3a) P-sep(3b) TS(3) HF/6-31G(d) 7.3 86.7 32.4 -15.1 59.0 HF/6-31++G(d,p) 6.3 62.9 10.7 -20.8 - B3LYP/6-31G(d) 10.2 93.0 49.4 - - B3LYP/6-31++G(d,p) 8.4 61.2 18.8 - - AM1 16.2 80.2 32.0 -21.0 34.9 The dipole moments (evaluated in center-of-mass coordinates) provide a rough indication of the charge distribution. Table S3 compares the AM1 values with the reference dipole moments from HF and B3LYP. Table S3. Dipole moments (D) for substrates and products. Species (1) (1a) (1b) TS(1) (2) (2a) (3) (3a) (3b) HF/6-31G(d) 1.90 8.76 5.43 4.86 11.05 16.02 5.70 3.24 8.83 HF/6-31++G(d,p) 1.93 9.12 5.91 4.90 11.64 16.48 6.13 3.31 9.76 B3LYP/6-31G(d) 1.79 7.73 4.03 4.24 10.24 12.99 4.91 3.21 - B3LYP/6-31++G(d,p) 1.97 8.34 4.91 4.28 11.12 14.19 5.35 3.34 - AM1 2.71 8.37 5.23 5.08 10.83 13.72 4.57 2.28 9.19 The AM1 geometries for selected reactants and products have already been compared to ab initio and B3LYP geometries in previous work [23] where good agreement has generally been found except for the O-O bond length. Table S4 presents additional comparisons for the geometry of the transition state in (1)+HOOH. The performance of AM1 is similar as before [23], the only major discrepancy is again the underestimation of the O-O distance (even though it should be noted that the increase of the O-O distance in going from HOOH to the transition state is well reflected). Table S4. Selected bond lengths (A) in the transition state for the reaction between the p-OHB dianion and hydrogen peroxide; see Figure 1 in the paper for the numbering of the atoms. Bond length C1-C2 C2-C3 C3-C4 C4-O4 C3-H3 C3-Od Od-Op Od-Hd Op-Hp HF/6-31G(d) 1.359 1.435 1.457 1.231 1.075 1.953 1.873 0.947 0.947 HF/6-31++G(d,p) 1.363 1.432 1.453 1.237 1.075 1.990 1.846 0.946 0.944 B3LYP/6-31G(d) 1.383 1.427 1.470 1.259 1.087 2.047 1.911 0.969 0.970 B3LYP/6-31++G(d,p) 1.390 1.420 1.461 1.269 1.085 2.209 1.815 0.970 0.968 AM1 1.380 1.415 1.447 1.267 1.098 1.965 1.646 0.960 0.950 In an overall comparison, AM1 reproduces the HF and B3LYP reference data in Tables S1-S4 with reasonable accuracy. References for this section: [S1] W.J. Hehre, L. Radom, P.v.R. Schleyer, and J.A. Pople, Ab Initio Molecular Orbital Theory, Wiley, New York, 1986. [S2] A.D. Becke, J. Chem. Phys. 98, 5648 (1993). [S3] Gaussian98, Revisions A.7 and A.9, M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, V.G. Zakrzewski, J.A. Montgomery, Jr., R.E. Stratmann, J.C. Burant, S. Dapprich, J.M. Millam, A.D. Daniels, K.N. Kudin, M.C. Strain, O. Farkas, J. Tomasi, V.Barone, M. Cossi, R. Cammi, B. Mennucci, C. Pomelli, C. Adamo, S.Clifford, J. Ochterski, G.A. Petersson, P.Y. Ayala, Q. Cui, K.Morokuma, D.K. Malick, A.D. Rabuck, K. Raghavachari, J.B. Foresman, J. Cioslowski, J.V. Ortiz, A.G. Baboul, B.B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R.Gomperts, R.L. Martin, D.J. Fox, T. Keith, M.A. Al-Laham, C.Y. Peng, A. Nanayakkara, C. Gonzalez, M. Challacombe, P.M. W. Gill, B. Johnson, W. Chen, M.W. Wong, J.L. Andres, C.Gonzalez, M. Head-Gordon, E.S. Replogle, and J.A. Pople, Gaussian, Inc., Pittsburgh PA, 1998. 2. Unrestrained environment in MD simulations ---------------------------------------------- The following parts of the system were unrestrained during the entire MD simulation. a) QM region (substrate and oxidised cofactor FADHOOH). b) 75 amino acids: 43 ILE, 44 ARG, 45 ALA, 46 GLY, 47 VAL, 48 LEU, 49 GLU, 70 LEU, 71 VAL, 72 HIS, 73 GLU, 74 GLY, 75 VAL, 76 GLU, 77 ILE, 78 ALA, 79 PHE, 82 GLN, 83 ARG, 84 ARG, 85 ARG, 86 ILE, 88 LEU, 89 LYS, 92 SER, 95 LYS, 96 THR, 97 VAL, 98 THR, 99 VAL, 101 GLY, 185 TRP, 187 GLY, 194 PRO, 196 SER, 197 HIS, 198 GLU, 199 LEU, 200 ILE, 201 TYR, 202 ALA, 203 ASN, 208 PHE, 209 ALA, 210 LEU, 211 CYS, 212 SER, 213 GLN, 214 ARG, 220 ARG, 221 TYR, 222 TYR, 247 LEU, 248 PRO, 292 PRO, 293 PRO, 294 THR, 295 GLY, 296 ALA, 297 LYS, 300 ASN, 342 PHE, 346 MET, 350 LEU, 378 LEU, 379 ALA, 380 THR, 381 ILE, 382 ALA, 383 GLU, 384 ASN, 385 TYR, 386 VAL, 387 GLY, 388 LEU. c) 19 water molecules: 22 WAT, 23 WAT, 25 WAT, 31 WAT, 33 WAT, 34 WAT, 40 WAT, 64 WAT, 75 WAT, 90 WAT, 93 WAT, 103 WAT, 107 WAT, 150 WAT, 177 WAT, 178 WAT, 179 WAT, 181 WAT, 219 WAT. All listed residues (species) have at least one atom at a distance of less than 7 A from one of the centres of mass of the imaginary chain formed by the species: p-OHB, 201 TYR, 385 TYR, 33 WAT, 219 WAT, and 72 HIS. The numbering scheme comes from reference 15. 3. Definition of the MM' interaction region in MD simulations ------------------------------------------------------------- The MM' interaction region contains during the entire MD simulation: a) 63 amino acids: 12 PRO, 42 ARG, 43 ILE, 44 ARG, 45 ALA, 46 GLY, 47 VAL, 48 LEU, 49 GLU, 52 MET, 72 HIS, 75 VAL, 97 VAL, 98 THR, 99 VAL, 100 TYR, 101 GLY, 102 GLN, 184 GLY, 185 TRP, 186 LEU, 187 GLY, 188 LEU, 189 LEU, 199 LEU, 200 ILE, 201 TYR, 208 PHE, 210 LEU, 211 CYS, 212 SER, 213 GLN, 214 ARG, 219 SER, 220 ARG, 221 TYR, 222 TYR, 262 GLU, 264 SER, 265 ILE, 266 ALA, 267 PRO, 268 LEU, 286 ASP, 289 HIS, 290 ILE, 291 VAL, 292 PRO, 293 PRO, 294 THR, 295 GLY, 296 ALA, 297 LYS, 298 GLY, 299 LEU, 300 ASN, 301 LEU, 302 ALA, 303 ALA, 343 SER, 346 MET, 347 THR, 385 TYR. b) 25 water molecules: 13 WAT, 22 WAT, 24 WAT, 31 WAT, 33 WAT, 49 WAT, 58 WAT, 64 WAT, 68 WAT, 70 WAT, 71 WAT, 90 WAT, 92 WAT, 93 WAT, 103 WAT, 107 WAT, 128 WAT, 172 WAT, 173 WAT, 178 WAT, 179 WAT, 181 WAT, 198 WAT, 206 WAT, 219 WAT. The numbering scheme comes from reference 15. 4. Technical details of QM/MM geometry optimisations in PHBH ------------------------------------------------------------- The protocol for the QM/MM geometry optimisations of the full enzyme consisted of four phases: - MM and QM/MM relaxation of the initially prepared structure, - microiterative transition state (TS) search at the QM/MM level, - QM/MM energy minimisation for the reactant complex starting from TS, - QM/MM energy minimisation for the product complex starting from TS. Each of the four phases required a number of separate steps which are outlined below. The original paper on the HDLC optimiser [38] should be consulted for further background information. Relaxation at the MM level (a-f) and the QM/MM level (g): a) Starting point: the initially prepared structure (see section III of the paper), partitioned into a central region (p-OHB, FADHOOH, Arg214, Tyr201, Tyr385) and the remainder ("environment"). b) Energy minimisation of the environment with fixed central region using 40 HDLCopt cycles. c) Energy minimisation of the central region with fixed environment using 20 HDLCopt cycles which defines the reference positions of the central region for the next steps. d) 5 ps MD simulation and subsequent steepest descent minimisation of the entire system with harmonic position restraints in the central region, using a restraining force constant of 100 kcal/(mol*A**2). e) Repeat previous step with a force constant of 50 kcal/(mol*A**2). f) Repeat previous step with a force constant of 25 kcal/(mol*A**2). g) Refine the resulting structure by an energy minimisation using the HDLC optimiser with the input option cfact=0.5 and a convergence criterion for the gradient components of gtol=0.002 Hartree/Bohr: this leads to a reasonably strain-free reactant complex since the starting point (a) was derived from the X-ray structure of a PHBH-substrate complex [15]. Transition state search at the QM/MM level: h) Define a reaction core of 9 atoms (see Figure 5 of ref.38). i) Shift distal OH group in the relaxed structure (g) by 0.1 A towards the C3 atom of the substrate. j) Energy minimisation of all HDLC fragments except for the reaction core using 100 HDLCopt cycles. k) Energy minimisation of the HDLC fragments in the QM region outside the reaction core for another 100 HDLCopt cycles. l) Microiterative transition state search in the QM region with gtol=0.00135 Hartree/Bohr followed by an energy minimisation of all HDLC fragments (except for the reaction core) with gtol=0.00045 Hartree/Bohr. m) Repeat previous step twice. n) Final microiterative transition state search with gtol=0.0009 Hartree/Bohr for the reaction core and gtol=0.0003 Hartree/Bohr for the other HDLC fragments. QM/MM energy minimisation of the product complex: o) Shift distal OH group from optimised TS geometry (n) by 0.1 A towards C3 atom of substrate. q) HDLCopt energy minimisation with gtol=0.00045 Hartree/Bohr. QM/MM energy minimisation of the reactant complex: r) Shift distal OH group from optimised TS geometry (n) by 0.1 A towards Op atom of cofactor. s) HDLCopt energy minimisation with gtol=0.00045 Hartree/Bohr. 5. Technical note concerning Figure 5 of the paper -------------------------------------------------- The originally generated profile [21] for monoanion 2 with model A (Figure 5) showed some structure in the product region, for values of the reaction coordinate R(Od-C3) below 1.4 A. We have checked the original evaluations again and found that there were unusually large variations in dG/dlambda during the course of the simulations for two of the corresponding lambda points. We now consider these data as being technically unreliable and have therefore not included them in Figure 5. This omission does not affect any of our conclusions.