Tuesday, March 15, 2011

Discussion of the Thiel paper

The discussion page for the Thiel paper:

'Coupling and uncoupling mechanisms in the methoxythreonine mutant of cytochrome P450cam: A quantum mechanical/molecular mechanical study.'

Collection of the various points mentioned. Suggestions, feedback, corrections, additions and comments are most welcome.

General notes on the paper:

  • The present mutagenesis study wishes to investigate the effect of the Thr252MeO-Thr mutation in terms of disruptions of the proton relay channels (Asp251 and Glu366) which are essential prerequisites for the conversion of Cpd0 into Cpd1 in the cytochrome P450cam enzymatic cycle. The article thus continues the line of computational mutagenesis research studies initiated by the same group in 2009 [1].
  • The study performs a single classical MD simulation from an initial PDB structure (1DZ8) obtained from group of Schlichting et al. [2] and a single (random) snapshot is taken as initial structure for the QM/MM optimizations. Hydrogen link atoms and the charge shift model were employed to treat the QM/MM boundary and an electronic embedding scheme [3],[4] was adopted in the QM/MM calculations. In the QM calculations, interactions with MM charges were incorporated into the LCAO-MO one-electron Hamiltonian, thus allowing QM polarization, and the QM/MM electrostatic interactions were evaluated from the QM electrostatic potential (at HF/6-31G(d) level) and the MM partial charges.
    The procedure is hence sequential; the QM region is initially optimized with the surroundings fixed (1 QM step) followed by a MM response to the QM region which is now being held fixed. This is again followed by a new QM step (energy and gradient evaluation for the core) and the iteration cycle continues until convergence [5]. Transitions states may thus be located from the inversion of the QM explicit Hessian by the P-RFO module of the HDCL procedure within ChemShell. P-RFO updates the Hessian matrix via the quasi-Newton BFGS optimization algorithm. The QM calculations were done with UB3LYP in Turbomole as a compromise between the need for incorporation of static electronic correlations effects and the large size of the QM region while all the MM calculations were done with the CHARMM22 force field either in CHARMM (MD) or in DL-POLY.
  • The results of the study show that the computed QM/MM barriers indicate unfavorable uncoupling in the case of the Thr252MeO-Thr mutant (unlike when X = Val, Ala, Gly [1]), whereas there are two energetically feasible proton transfer pathways for coupling; these are Mechanism I: homolytic O-O bond cleavage followed by coupled proton-electron transfer and Mechanism II: proton-assisted heterolytic O-O bond cleavage. The study shows Mechanism I to be favorable.
    The corresponding rate-limiting barriers for the formation of Cpd1 are higher in the mutant than in the wild-type enzyme. These findings are consistent with the experimental observations that the Thr252MeO-Thr mutant forms the alcohol product exclusively (via Cpd1), but at lower reaction rates compared with the wild-type enzyme.
    With respect to the two different channels, the rate-limiting barriers are somewhat lower in the Glu366 channel than in the Asp251 channel. The Asp251 channel is in contact with bulk water, though, so it should be rather facile to reprotonate Asp251 after each coupling reaction that involves proton transfer in the Asp251 channel. This is not true for Glu366, which resides in a hydrophobic pocket and is thus difficult to reprotonate.

Topics of discussion:

  • We discussed the potential problems associated with only employing a limited number of MD simulation snapshots in the subsequent QM/MM calculations. This procedure surely does not incorporate effects affiliated with the dynamics of the system, i.e. effects arising from conformational fluctuations are neglected as the span of the configurational phase space is largely limited. The Thiel group includes a short assessment of this issue as they address another (randomly chosen) snapshot in the same manner and achieve similar results, thus confirming the validity of the procedure.
  • At the lecture we discussed this aproval of validity of the computational protocol - which is also build upon the excellent results with respect to experimental references which the group has obtained in previous studies [6],[7] - and we all felt that the neglect of dynamics in terms of the limited amount of snapshot used may potentially cause problems but we all acknowledge the apparent accuracy of the method. We spoke further on the neglect of solvent effects as the hydrogen bonding pattern in the solvent might change dramatically within the time period of the MD simulation potentially leading to increments in energy barriers. In defence of the present study, the active site of cytochrome P450cam resides within a hydrophobic pocket of the enzyme but fluctuations (+/-) in the energy barriers as a result of the shifting hydrogen bonding patterns in the solvent may still occur.
  • Finally, we questioned wether the Pople style basis set, 6-31G, i.e. without polarizations functions, would be efficient in describing the electronic environment within the QM calculations.

Selected references:

[1] Altarsha et al. - J. Am. Chem. Soc., 2009, 131, 4755
[2] Schlichting et al. - Science, 2000, 287, 1615
[3] Altun et al. - J. Phys. Chem. B, 2005, 109, 1268
[4] Bakowies et al. - J. Phys. Chem., 1996, 100, 10580
[5] Turner et al. - Phys. Chem. Chem. Phys., 1999, 1, 1323
[6] Schöneboom et al. - J. Am. Chem. Soc., 2002, 124, 8142
[7] Schöneboom et al. - J. Phys. Chem. B, 2004, 108, 7468

For further reading, please see the list of references at the end of the slide show.


Jan Jensen said...

I didn't see a link to the original paper.

Another thing we discussed was whether the study uses the conventional B3LYP functional (as the referencing in the paper suggests) or the reparameterized version you mention on slide 10.

Also, which of your references describe the solvation methodology you mention on slide 8?

Finally, the discussion is part of a graduate course on QM/MM.

Janus J. Eriksen said...

@Jan >

a) Yes, sorry for the missing links and thanks for the rectifications.

b) Wrt the solvation methodology (protonation states, etc.), this is most detailed described in Ref. [7], also refered to as Schöneboom2004 in the presentation.

c) Yes, we discussed the reasoning behind the choice of DFT functional. The group uses a unrestricted KS DFT formalism as the 'correct' state of the Fe-O compound (Cpd0) is in fact a spin contaminated state as it is a mixture of a doublet state (the most important in an energetic sense) and a quartet and sextet state (see [6]/Schöneboom2002 and [7]). There is thus a need for being able to describe the static electron correlation effects and UB3LYP (with the original 20% HF B3 exchange) is chosen as application of WFT MSSCF is not realistic on a molecular system of this size. Ref. [7] compares B3 with B88 (no HF exchange) and the Becke half-and-half functional (50% HF exchange) and Ref. [3] and the supporting material of [3] discusses the importance of HF admixture in detail.