Does Proton Conduction in the Voltage-Gated H+ Channel hHv1 Involve Grotthuss-Like Hopping via Acidic Residues?

Siri C. van Keulen, Eleonora Gianti, Vincenzo Carnevale, Michael L. Klein, Ursula Rothlisberger, and Lucie Delemotte (2017)
Contributed by Dries Van Rompaey

Voltage gated proton channels are membrane channels that are are regulated by both the pH gradient and the voltage. In most cases these channels only open when the electrochemical proton gradient is outwards, functioning as a passive acid extrusion mechanism. They have an exquisite selectivity for protons, prompting Delemotte and coworkers to investigate the mechanism of proton transport and the origin of this selectivity.

The mechanism for proton transport was explored using an integrative approach, combining classical simulations with QM/MM. Three separate cation binding sites can be observed in the channel, each consisting of a pair of negatively charged residues. QM/MM simulations were used to examine binding and unbinding of the proton, while conformational rearrangements were sampled using classical simulation. In contrast to the classic Grotthus mechanism, where the proton is transported through a water wire by subsequent bond formation and bond breaking, simulations indicated that the proton may traverse the channel by jumping between the acidic residues, assisted by a water molecule. As shown by classical MD, structural rearrangements occur upon proton binding, orienting the residues for the next proton jump, allowing for the proton to move through the channel. This integrative modelling study provides an interesting and novel mechanism for the transfer of a proton across a membrane. It will be interesting to see where future work on this channel leads.

Multiscale Quantum Mechanics/Molecular Mechanics Simulations with Neural Networks

Lin Shen, Jingheng Wu, and Weitao Yang (2016)
Contributed by Jan Jensen

There has been a lot of work on estimating high level energies from low level energies, most of which have focussed some kind of interpolation or extrapolation.  However, this has been particularly challenging for QM/MM-MD studies where thousands of semiempirical energies contribute to the PMF but only hundreds of high level calculations are feasible. 

Yang and co-workers offer an interesting solution to this problem by using the high level calculations to train a neutral network to estimate the energy correction, which can then be used to estimate the high level energies for all thousands of geometries from the semiempirical QM/MM-MD.  Thus, the PMF can be computed "from scratch" at the high level rather than correcting the low level PMF.

The approach is based on work by Behler and Parrinello, but with three tweaks geared towards this particular problem.

1. The neural network is trained to reproduce an energy correction rather than a total energy.

2. Mulliken charges are used as input, in addition to atomic coordinates of the QM region, to provide some information about the interaction with the MM environment.

3. A subnet is added to the neutral net that depends on the reaction coordinate and the potential of mean force obtained at the low level to improve accuracy.  

The method is tested for three relatively small systems in water, so that high level PMFs can be computed rigorously for comparison. The results are quite impressive. For example, the free energy of glycine zwiterion is ca 8 kcal/mol higher in energy than the neutral form at the SCC-DFT/MM level but ca 8 kcal/mol lower in free energy at the B3LYP/6-31G(d)//MM level of theory. This latter value is reproduced to within 0.1 kcal/mol by the machine learning approach using only 500 B3LYP/6-31G(d)//MM energies to train the neural net. For comparison, the PMF requires 50,000 energy evaluations.

This work is licensed under a Creative Commons Attribution 4.0

Chromophore−Protein Coupling beyond Nonpolarizable Models: Understanding Absorption in Green Fluorescent Protein

Csaba Daday, Carles Curutchet, Adalgisa Sinicropi, Benedetta Mennucci, and Claudia Filippi, 
J. Chem. Theory Comput. 2015, 11, 4825 – 4839
Contributed by Tobias Schwabe

Modeling enzymichromism, the sprectral tuning of a chromophore by its protein surroundings, is a formidable challenge for computational chemists. It involves methods from molecular dynamics simulations to excited state multi-reference computations in a QM/MM framework (or even more advanced embedding schemes).

Obviously, in each step of the modeling one has to choose from several methods and different applications. As a consequence, no standard protocol how to address the problem has been agreed upon so far. Daday et al. now presented a thorough analysis of possible routes and identified key ingredients for success in enzymichromism modeling in the paper highlighted here. As a test case, they picked the green fluorescent protein (GFP) and the electronic excitation of its chromophore in the protonated (A) and deprotonated (B) form. Two main issues are covered: What is the effect of using an optimized protein structure vs. snapshots from a MD run (and average) and how should the protein environment be modeled?

Let's look at the MD results first: For an average of 50 snapshots, the excitation energies are 3.36 ± 0.1 eV (A form) and 3.07 ± 0.1 eV (B form) at the CAM-B3LYP/6-31+G*/MM level of theory, but the energy spread is 0.5 eV in both cases. Nevertheless, the average is close to the value of an annealed structure (representative for an optimized crystal structure), which the authors also computed: 3.38 eV (A form), 3.11 eV (B form). It is hard to tell if this finding might be transferable to other protein systems or if this good agreement between ensemble average and optimized structure is just a lucky match. The energy spread might be a caveat.

In any case, the obtained values are blue-shifted in comparison to experiment which is not a failure of CAM-B3LYP (alone). This has been checked by recalculating the results with CASPT2/MM. Therefore, the authors also tried different embedding schemes instead of static point charges: state-specific induced dipoles, linear response methods including induced dipoles, and even frozen density embedding (Those readers who want to learn more about the methods and their differences, should refer to [1]). The bottom line: frozen density embedding is not a significant improvement over state-specific induced dipoles, but combining the latter with the effects of linear response in a polarizable embedding yields very good results in comparison. This is in accordance with our previous study of the problem for which CC2 in a polarizable embedding (PERI-CC2) has been applied.[2] There, a very good agreement between PERI-CC2 and all-QM computations has been demonstrated (for smaller cluster models, of course).

The study highlighted here emphasizes the importance of coupling induced dipoles to the transition moments of the excitation which is not done in all methods which allow for polarization in the classical region of a QM/MM approach. These findings also concern the even broader field of solvation models like PCM or COSMO. Everyone interested in environmental effects on electronic excitations (and other dynamical properties) should make oneself familiar with the differences of state-specific and linear response treated polarization effects. Often, they are not pointed out clearly in the literature which hampers a good interpretation of computed results and their assessment in comparison to experiment.


[1] A. S. P. Gomes, C. R. Jacob, Quantum-chemical embedding methods for treating local electronic excitations in complex chemical systems, Annu. Rep. Prog. Chem.,Scet. C:Phys. Chem. 2012, 108, 222–277 (http://dx.doi.org/10.1039/c2pc90007f)
[2] T. Schwabe, M. T. P. Beerepoot, M. T. P., J. M. H. Olsen, J. Kongsted, Analysis of computational models for an accurate study of electronic excitations in GFP Phys. Chem. Chem. Phys. 2015, 17, 2582–2588. (http://dx.doi.org/10.1039/C4CP04524F)

Computational Design of a Time-Dependent Histone Deacetylase 2 Selective Inhibitor

Jingwei Zhou, Min Li, Nanhao Chen, Shenglong Wang, Hai-Bin Luo, Yingkai Zhang , and Ruibo Wu

Histone deacetylases (HDACs) are a family of enzymes involved in gene expression and post-translational modifications. HDACs are very important targets for drug development due to their roles in cancer and other diseases. Several HDAC inhibitors have been developed, however, many known inhibitors produce side-effects because of their poor selectivity. This paper presents a striking example of how state-of-the-art QM/MM-MD calculations can be used for computer-aided drug design (CADD). 

Based on their previous calculations on the mechanism of the wild-type enzyme, Zhou and co-workers hypothesized that a selective inhibitor could be created by developing a molecule that would undergo an HDAC-catalyzed intra-molecular reaction. To this end, the authors proposed several candidates and used QM/MM-MD simulations to study the mechanism in gas-phase, solution and in the enzyme active site. The most promising candidates ( –hydroxymethyl and –aminomethyl substituted chalcones) were subsequently synthesized and characterized in vivo. The authors further analyzed the mechanism of inhibition via QM/MM-MD for the most promising candidate to understand how this molecule acts as an HDAC2-selective, time-dependent inhibitor.

 Reprinted with permission from ACS Chem. Biol., Article ASAP DOI: 10.1021/cb500767c . Copyright (2014) American Chemical Society.

Comparison of Ab Initio, DFT and Semiempirical QM/MM approaches for description of catalytic mechanism of hairpin ribozymes

Vojte Mlynsky, Pavel Banas, Jiri Sponer, MarcW. van der Kamp, Adrian J. Mulholland, and Michal Otyepka, Journal of Computational and Theoretical Chemistry, 2014, DOI:10.1021/ct401015e

Contributed by Esteban Vohringer-Martinez

I enjoyed very much reading this paper where the authors report a very interesting study in which the performance of semiempirical methods in QM/MM calculations is addressed in detail.

The reaction under study is the reversible phosphodiester bond cleavage and ligation in hairpin ribozymes. These ribozymes belong to a small group which catalyze the reaction without any metal ion at comparable rates.

From experimental evidence and previous theoretical calculations it was known that there are two main players in the catalytic reaction: A38 and G8; two neighbouring bases to the phosphodiester bond to be cleaved between the Adenine -1 and Guanine +1 (see scheme below). However which was not clear from the theoretical studies nor the experiment is the protonation state for the A38 and G8 bases. The experimental studies showed a pH dependence in the catalytic activity of the ribozyme implying different catalytic activity as a function of protonation state.

To account for different protonation states the authors studied two reaction mechanisms: the mono anionic one (top) where only the phosphatediester group is deprotonated and the dianionic mechanism where the G8 bases is also deprotonated (bottom). In both mechanism they assumed the A38 to be protonated (see reaction scheme).

The energetics of the reaction are followed along two reaction coordinates represented by a linear combination of the cleaving and forming bonds (d1 and d2) in the proton transfer step and the oxygen-phosphor distance as shown in the scheme. 

The energies were calculated with the SCS(spin component scaled)-MP2  single point energy calculations (CBS limit) on BLYP optimized geometries (6-31G(d,p) in a QM/MM electrostatic embedding with the AMBER99sb force field.  The test methods include the MPW1K hybrid functional, the BLYP GGA DFT functional and the semiempirical methods AM1/dPhoT and SCC-DFTBPR.

For me the main finding of this paper is that Ab-Initio and all DFT methods predict a concerted nucleophilic and proton transfer mechanism whereas both semiempirical methods yield a sequential mechanism where first the proton is transferred and then the nucleophilic attack takes place. Interestingly both semiempirical methods yield a stable intermediate in the dianionic mechanism after the proton transfer step which is not present in the ab-initio and DFT results. 

These result raise doubts about the ability of semiempirical methods to provide correct structures and reaction mechanisms in enzyme catalysis and suggest that a careful calibration of these methods for each system should be performed prior to its usage. 

One additional aspect to consider is that the activation barrier of the semiempirical methods were close to the ab-initio, DFT methods and experiment. However, a match of barrier heights as shown in this study does not guarantee a correct reaction mechanism.

Finally, the authors also compare their results to free energy calculations employing semiempirical methods. But, as the authors also conclude the almost negligible entropy contribution they report should be taken with care due to the very short sampling time of only 50ps.

Computer Simulation and Analysis of the Reaction Pathway of Triosephosphate Isomerase

P. A. Bash, M. J. Field, R. C. Davenport, G. A. Petsko, D. Ringe, and M. Karplus Biochemistry 1991, 30, 5826

Contributed by +Jan Jensen 

In honor of the 2013 Nobel Prize in Chemistry I am highlighting this gem from 22 years ago by Karplus and co-workers.  I believe (correct me if I'm wrong) that this paper presents the first study of enzyme catalysis using "conventional" QM/MM: conventional electronic structure theory (AM1) combined with a standard protein force field (CHARMM).    

When the paper came out I had recently started as a PhD student and one of my projects was working on the Effective Fragment Potential QM/MM method.  At that time it was relatively easy to keep track of the QM/MM literature and two papers were usually on top of my rather short stack of QM/MM papers: Singh & Kollman and Field, Bash & Karplus.  As I remember it, I was reading them during an extended visit to the Center for Advanced Research in Biotechnology in Maryland and sitting there it was pretty hard to imagine, based on the applications in these papers, how QM/MM would ever help advance research in biotechnology.  There was a hint in the Field, Bash, and Karplus paper, where they used something called Triose Phosphate Isomerase (TIM) to motivate the use of link atoms, but it wasn't until the Biochemistry paper that all the pieces were put together.  

All of a sudden (semi-empirical) electronic structure theory could be used to say something meaningful about a system with a thousand atoms (1650 to be exact).  Specifically that Lys12 is most important to catalysis and His95 could act as an acid/base catalyst.

Despite being the first of its kind, the paper still reads very much like a "modern" study and I think it is fair to say that most subsequent QM/MM studies of enzymes are really variations on the themes introduced in this paper.

This work is licensed under a Creative Commons Attribution 3.0 Unported License.

Unraveling the Enigmatic Mechanism of L‐Asparaginase II with QM/ QM Calculations

Diana S. Gesto, Nuno M. F. S. A. Cerqueira, Pedro A. Fernandes, and Maria J. Ramos J. Am. Chem. Soc. 2013, 135, 7146

Contributed by Jonathan Goodman

This paper discusses the hydrolysis of asparagine

The transformation is simple, but the catalyst, L-Asparaginase II, is not. The active site of the enzyme contains a number of functional groups which work together to hydrolyse the enzyme. 

The currently accepted mechanism is that the process is related to that of a serine protease, with Thr12 attacking the amide to form an ester that is then hydrolysed. This mechanism looks reasonable. The three-dimensional structure shows that Thr12 is poised over the carbonyl of the amide, close to the π* orbital of the carbonyl and in about the right position to attack. However, studies of mutants of the enzyme do not provide unambiguous support for this proposal. 

This paper reports a computational study of the system using an ONIOM approach, with B3LYP/6-31G(d) for the high-level layer and AM1 for the low-level layer. This split approach made it possible to follow the reaction pathways for the system. Minima and transition states were recalculated using single point calculations at the M06-2X/6-311++G(2d,2p) level. This method showed that the currently accepted mechanism, summarized in the figure below showing the key Thr12-substrate interaction in cyan, was a high energy pathway, and an alternative mechanism involving the attack of a water molecule and stabilization of the intermediate oxyanion by Thr12 was more accessible. The new mechanism is surprising in that Thr12 is above the carbonyl group and so not in a good position to interact with the lone pair of the amide carbonyl. However, this configuration does fit the grand jeté orientation that we recently highlighted for mechanisms involving oxyanion holes (doi: 10.1039/C2OB06717J).

This paper demonstrates that the calculations can provide a sufficiently accurate analysis of the system to inspire an alternative to the currently accepted mechanism. The new mechanism fits all of the available experimental data; it was hard to see how some of the experimental data fitted the original mechanism. The new mechanism, therefore, is to be preferred. If both mechanisms had fitted the available experimental data, would these demanding calculations have been enough to change the currently accepted mechanism to the new process?  

Dispersion corrections and bio-molecular structure and reactivity

Richard Lonsdale, Jeremy N. Harvey, and Adrian J. Mulholland "Effects of Dispersion in Density Functional Based QM/ MM Calculations on Cytochrome P450 Catalysed Reactions"Journal of Chemical Theory and Computation 2012, ASAP (Paywall)

Takeshi Nagata, Dmitri G. Fedorov, Hui Li, and Kazuo Kitaura "Analytic gradient for second order Møller-Plesset perturbation theory with the polarizable continuum model based on the fragment molecular orbital method" Journal of Chemical Physics 2012, 136, 204112 (Paywall)

The DFT dispersion correction developed by Grimme and co-workers was recently highlighted in Computation Chemistry Highlights. Recently two papers have appeared that quantify the importance of the dispersion correction on modeling bio-molecular structure and reactivity.

Cytochrome P450 barrier heights in better agreement with experiment using B3LYP-D2 and -D3

Lonsdale et al. used QM/MM to compute barrier heights for oxidation reactions, catalyzed by P450$_{cam}$, with an without dispersion corrections in the QM region.  Invariably the dispersion correction lowered the barrier significantly (usually by ca 5 kcal/mol), yielding results that were in better agreement with experimental values.  The effect of the dispersion correction on the transition state geometries was less pronounces though not negligible, with bond lengths changing by as much as 0.2 Å.

It is worth noting that the QM region contains a conjugated porphyrin ring and that three of the four substrates considered in the study contain one or more double bonds.  Thus, the QM region contains very polarizable functional groups and, since the magnitude of dispersion interactions increases with the polarizability of the groups involved, it is possible that the effect of dispersion on barrier heights for other enzymes will be less than observed here.  It will be interesting to find out.

MP2 quality Trp-cage structure using RHF-D

Nagata et al. have implemented analytical MP2/PCM gradients for the fragment molecular orbital method and used it to geometry optimize the small protein Trp-cage at the MP2/6-31(+)G(d) level of theory [the (+) indicates diffuse functions on carboxylate groups].  The resulting structure compared well with the experimental NMR structures, with a backbone RMSD of only 0.426 Å. This is a significant improvement in agreement compared to the corresponding RHF/PCM optimized structure (RMSD 1.107 Å) and demonstrates the importance dispersion in bio-molecular structure.  Interestingly, the corresponding RHF-D structure compared equally well to experiment (0.414 Å) and was virtually identical to the MP2 structure (RMSD 0.068 Å).

Disclaimer: I was involved in the implementation of the FMO RHF/PCM interface.

This work is licensed under a Creative Commons Attribution 3.0 Unported License.