A Local Variant of the Conductor-Like Screening Model for Fragment-Based Electronic-Structure Methods

Albrecht Goez and Johannes Neugebauer J. Chem. Theory Comput., Article ASAP, DOI: 10.1021/acs.jctc.5b00832
Contributed by Christoph Jacob
 
Fragment-based methods nowadays make it possible to perform quantum-chemical calculations for rather large biomolecules, for instance light-harvesting protein systems [1]. Such methods are based on the idea of splitting a protein into smaller fragments, such as its constituting amino acids. This leads to a linear scaling of the computational effort with the size of protein. Popular examples of fragment-based electronic structure methods include the fragment molecular orbital (FMO) method [2] and a generalization of the frozen-density embedding scheme (3-FDE) [3].

In a recent article in JCTC, Goez and Neugebauer from the University of Münster (Germany) address an additional bottleneck that appears in such calculations. Usually, it is necessary to include a solvent environment in the calculations, in particular if charged amino-acid side chains are present. The simplest way of doing so are continuum solvation models, such as COSMO or PCM. These models represent the solvent in terms of apparent charges on the surface of a cavity enclosing the protein. However, for proteins the number of apparent surfaces charges becomes rather large - for ubiquitin, a protein with only 78 amino acids, already 20,000 charges are needed. Updating these apparent surface charges involves solving a linear system of equations of size 20,000 x 20,000. When doing so in each SCF cycle for each of the fragments, the continuum solvation model will become the bottleneck of the calculation.

To solve this problem, Goez and Neugebauer developed a local variant of the COSMO model (LocCOSMO). In each fragment calculation, they update only those apparent surface charges that are close to this fragment. This reduces the computational effort significantly, but because every fragment is updated at some point it will eventually result in the same final result. This is demonstrated by the authors for several test cases. They can reduce the computational time required for a 3-FDE calculation of ubiquitin in a solvent environment by a factor of 30, without compromising the quality of the result.

The efficient combination of fragment-based quantum chemistry with continuum solvation models  provides an important tool for studies of biomolecules. It will make such calculations more robust by alleviating convergence problems for charged amino acids and will allow for a more realistic inclusion of protein environments in studies of spectroscopic properties of chromophores in biomolecular systems.

References:

[1] A. Goez, Ch. R. Jacob, J. Neugebauer, “Modeling environment effects on pigment site energies: Frozen density embedding with fully quantum-chemical protein densities”, Comput. Theor. Chem. 1040–1041, 347–359 (2014).
[2] D. G. Fedorov, K. Kitaura, “Extending the Power of Quantum Chemistry to Large Systems with the Fragment Molecular Orbital Method”, J. Phys. Chem. A 111, 6904–6914 (2007).
[3] Ch. R. Jacob, L. Visscher, “A subsystem density-functional theory approach for the quantum chemical treatment of proteins”, J. Chem. Phys. 128, 155102 (2008).

Simulations of Chemical Reactions with the Frozen Domain Formulation of the Fragment Molecular Orbital Method

Hiroya Nakata, Dmitri G. Fedorov, Takeshi Nagata, Kazuo Kitaura, and Shinichiro Nakamura. Journal of Chemical Theory and Computation 2015, 11, 3053−3064
Contributed by +Jan Jensen

Reprinted (adapted) with permission from J. Chem. Theory Comput. 2015, 11, 3053−3064. Copyright (2015) American Chemical Society

I am highlighting this paper because it is the first paper I have come across that reports the optimized TS, free energy barrier, and IRC for an enzymatic reaction using an all-QM description of the entire enzyme.

The underlying method is the frozen domain with dimers (FDD) formulation of the fragment molecular orbital method (FMO) implemented in the GAMES program.  The FMO is fragmentation method and in the FDD formulation the system is divided into three domain: A) an active domain where the both the geometry and density are allowed to change, B) a polarizable buffer where only the density is allowed to change and F) a frozen domain where neither is allowed to change.  The latter restraints means that the fragments in the frozen domain only have to be calculated once.

This paper presents the implementation of the fully analytic gradient and Hessian at the RHF level and the interface to the IRC methodology in GAMESS.  (The interface to the TS-optimization code was reported previously using a different variant of FMO.)  In this approach only the Hessian of the active domain is computed.

The method is then applied to the keo-enol tautomerization of phosphoglycolohydroxamic acid (PGH) by triosephosphate isomerase (TIM).  The active domain is PGH and the proton donor Glu410 while the polarizable buffer is any residue or water molecule within 5.2 Å of the active domain, which were both treated at the RHF/6-31G(d) + dispersion level of theory.  The frozen domain, i.e. the rest of the TIM monomer, is treated at the RHF/STO-3G level of theory and the entire system is comprised of 9227 atoms.

The workflow was at follows using 128 2.93 GHz Xeon cores:

1. Constrained optimization to obtain TS guess (12 hrs)
2. Hessian calculation (24 hrs)
3. TS optimization (12 hrs)
4. Hessian calculation (24 hrs)
5. IRC calculations (16 points) in each direction (42 hrs)
6. Optimization of reactants and products starting from last IRC points (12 hrs each)
7. Hessian calculation for reactant and products (24 hrs each)

Of course, at this level of theory and with such a small active domain the resulting free energy barrier (21.4 kcal/mol) is considerably different from the experimental value (14 kcal/mol) but this is only a proof-of-concept study.

An important next step for FMO studies of enzyme catalysis will be the implementation of a mixed DFT/RHF-D model where all or part of the active domain can be treated with DFT.  The computational cost-increase associated with using DFT could be off-set by using, for example, HF-3c for the non-reactive part of of the active region as well as the buffer region since this method gives excellent structures as low computational cost.

I thank +Casper Steinmann for alerting me to this article


This work is licensed under a Creative Commons Attribution 4.0

Three- and four-body corrected fragment molecular orbital calculations with a novel subdividing fragmentation method applicable to structure-based drug design

Chiduru Watanabe, Kaori Fukuzawa, Yoshio Okiyama, Takayuki Tsukamoto, Akifumi Kato, Shigenori Tanaka, Yuji Mochizuki, and Tatsuya Nakano Journal of Molecular Graphics and Modelling 2013, 41, 31.
Contributed by +Jan Jensen


One of the advantages of fragmentation based electronic structure methods over conventional linear scaling method is that one gets an automatic energy analysis in terms of interfragment interaction energies (IFIEs).  When applied to protein-ligans interactions such IFIE analyses can aid in the drug design process.  

However, in order to obtain accurate results the whole ligand is often treated as one fragment and the contribution of individual ligand functional groups to the binding is lost.  A main reason that accuracy depends on fragment size is that as fragments get smaller, many-body short-range interaction energies such as exchange repulsion become important.  Most fragmentation methods only treat such effects at the two-body level, while the Fragment Molecular Orbital (FMO) method also has a three-body corrected version (FMO3).  Nakano et al. have recently developed a four-body version of the FMO4 which yields accurate total energies for functional group-sized fragments, and this study by Watanabe et al presents the corresponding IFIE analysis.

At the MP2/6-31G level of theory and using an FMO2 calculation using one residue per monomer as a reference, the accuracy and CPU cost is comparable to an FMO3 calculation if the fragment size is roughly halved.  The CPU increases only by a factor of three on going to FMO4, while the accuracy is increased by an order of magnitude.

The relatively modest increase in CPU time is due to the fact that the many-body short-range corrections only are computed for fragments in close proximity to one another. Thus, the number of corresponding trimer and tetramer calculations scale roughly linearly with system size.  I have explored the question of computational efficiency and fragment-size further in this post.


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

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.  

Energy Decomposition Analysis in Solution Based on the Fragment Molecular Orbital Method

D. G. Fedorov, K. Kitaura, J. Phys. Chem. A 116 (2012) 704-719.


The solvent screening has long been described relying on the dielectric constants taken arbitrarily based on "experience". In this work the dielectric constants are derived from ab initio calculations using the charges induced by the solute on the solvent surface in the polarizable continuum model. These charges uniquely define the solvent screening energy, and also the dielectric constant.

Induced solvent charges δ surrounding two subsystems I and J.

[Reprinted with permission from the Journal of Physical Chemistry A.  Copyright (2012) American Chemical Society.]

For the rather trivial case of two interacting systems (sodium and chloride ions), the analysis provides rich physical insight into the solvation process, and the induced charge quenching effect is observed, caused by the partial loss of the solvation shell and the potential pressure of the nearby ion.
The real challenge for a physical interpretation is found in systems with multiple fragments, in which case strong many-body electrostatic effects lead to  very interesting and peculiar results, such as a partial solvent screening reduced by the pressure of the solute-induced electrostatic potential of nearby fragments.

The solvent screening of the electrostatic interactions between charged residues in a small protein chignolin (PDB: 1UAO) is discussed based on ab initio (FMO) calculations in solution.

Prediction of cyclin-dependent kinase 2 inhibitor potency using the fragment molecular orbital method

M. P. Mazanetz, O. Ichihara, R. J. Law, M. Whittaker Journal of Cheminformatics 2011, 3:2 (Open Access)

The paper describes a surprisingly accurate method for predicting ligand-protein binding energies based, in part, on MP2/6-31G(d) ligand-protein interactions computed using the Fragment Molecular Orbital (FMO) method.  The correlation between predicted and observed binding free energies are $r^2=0.939$ and $r^2=0.842$ for the training and testing set, respectively.  The method is sufficiently efficient to be used within an industrial setting.  The method performs better than more conventional scoring function implemented in the MOE program.

The binding free energy is given by$$\Delta G_{bind}=X_H \Delta H_{bind}+X_{TS} (T \Delta S_{bind})+X_{psolv} \Delta\Delta G_{psolv}+X_{npsolv} \Delta\Delta G_{npsolv}+X_{residual}$$
$\Delta H_{bind}$ is the ligand-protein interaction energy computed from an FMO energy decomposition analysis using a small structural model of the protein. $T\Delta S_{bind}$ is set to 1 kcal/mol per rotateable bond, $\Delta\Delta G_{psolv}$ and $\Delta\Delta G_{npsolv}$ are computed using the Poisson-Boltzmann equation and the solvent accessible surface area, respectively.

The four coefficients are obtained using a partial least squares to latent structures analysis on 14 ligand-protein complexes (for which x-ray coordinates were available).  Surprisingly, $X_{TS}$ came out negative, so it is doubtful whether this term actually reports on an entropy change. The method is then applied to a test set of 14 new ligands for which x-ray structures where not available.

Another interesting feature of the method is that only the hydrogen atoms of the x-ray structures were geometry optimized (using MMFF94x).  In the case where x-ray structures were not available, the structures of the ligand-protein complexes "were built by manually modifying the reference ligand and then energy minimising the ligand whilst keeping the reference heavy atoms fixed".

Results where  $\Delta H_{bind}$ is computed using the MMFF94x force field instead of FMO was not reported.  This would answer the question whether FMO is truly needed.

Acknowledgments: I thank Dmitri Fedorov for alerting me to this paper, and Michael Mazanetz for helpful discussions.  Note that in Table 2 the column labels for the experimental and computational binding free energies should be switched.

Update (4/3/2012): $X_H=0.018613, X_{TS}=-0.485403, X_{psolv}=-0.024463, X_{npsolv}=1.729698, X_{residual}=0.258$.  The latter parameter is in kcal/mol.

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