Quantum Chemical Approach to Estimating the Thermodynamics of Metabolic Reactions

Adrian Jinich, Dmitrij Rappoport, Ian Dunn, Benjamin Sanchez-Lengeling, Roberto Olivares-Amaya, Elad Noor, Arren Bar Even & Alán Aspuru-Guzik Sci. Rep. 2014, 4, 7022

Contributed by +Jan Jensen 

'

DOI: 10.1038/srep07022 CC-BY-NC-ND

Computational chemistry is finally in a position to conduct studies on hundreds or even thousands of molecules within a reasonable amount of time.  This has led to benchmark datasets which have been used to develop new semi-empirical methods (e.g. DFT-Dx or HF-3c) as well as high-throughput screening of molecules with desirable properties. 

This paper represents a third kind of study: a concerted effort to fill in missing experimental data needed to quantitatively understand a complex chemical system - in this case cellular metabolism. Standard Gibbs reaction free energies of biochemical reactions are crucial to quantitative modeling of metabolism (metabolic flux analysis) but only a small fraction have been measured experimentally.  

Many of the reactions involve molecule-sizes that can be modeled routinely by electronic structure calculations.  However, there are new challenges compared to, for example, atmospheric chemistry such as the modeling of pH- and solvation-effects and multiply charged ions.  Thus, any conformational search must often be repeated for several different protonation states and multiply charged ions often necessitate the treatment of explicit solvent molecules and, even then, result in relatively solvation energies with relatively large errors.  Also, multiply charged anions may present additional problems for DFT due to incomplete cancellation of the inter-electronic Coulomb repulsion at larger distances. 

This paper represents the first, pioneering, attempt to address these issues with high-throughput in mind and the results are encouraging. For reactions involving non-multiply charged species B3LYP/6-31G* with the COSMO solvation model augmented by explicit water molecules result in MADs less than 5 kcal/mol, while being fast enough to compute more than a hundred reaction energies (involving more than 5000 quantum calculations). 

However, this paper is more a call to action than the definitive answer.  Many new methodological improvements (perhaps such as this one?) are needed in order to increase the accuracy and efficiency needed produce a complete and accurate picture of cellular metabolism. I hope this paper will spur such developments - hopefully in the form of large and open collaborations (perhaps similar to this one?).  

Thanks to Adrian Jinich for alerting me to this paper

This work is licensed under a Creative Commons Attribution 4.0 International License.

Isotope Effects, Dynamic Matching, and Solvent Dynamics in a Wittig Reaction. Betaines as Bypassed Intermediates

Chen, Z.; Nieves-Quinones, Y.; Waas, J. R.; Singleton, D. A.  J. Am. Chem. Soc. 2014, 136, 13122-13125
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

If you hadn’t noticed, I am a big fan of the work that Dan Singleton is doing concerning the role of dynamics in discerning reaction mechanisms. Dan’s group has reported another outstanding study combining experiments, traditional QM computations, and molecular dynamics – this time on the Wittig reaction.¹

The key question concerning the mechanism is whether a betaine intermediate is accessed along the reaction (path A) or whether the reaction proceeds in a concerted manner (path B). Earlier computations had supported the concerted pathway (B).

Experimental determination of the heavy atom kinetic isotope effect was made for Reaction 1.

Reaction 1

Using the 6-31+G(2df,p) basis set, three different density functionals predict three different potential energy surfaces. With M06-2x, the surface indicates path A (stepwise), with the first step rate-limiting. B3P86 also predicts the stepwise reaction, but the second step is rate-limiting. The Lc-wPBE functional predicts a concerted reaction. Using these surfaces, they predicted the carbon isotope effect and compared it to the experimental values. The best agreement is with the M06-2x surface with a weighting of the vibrational energies of the two different TSs. The optimized structures of the two transition states, the betaine intermediate, and the product are shown in Figure 1.

TS1	Betaine
TS2	Product

Figure 1. M06-2x/6-31+G(2df,p) optimized geometry of the critical points of Reaction 1.

The agreement of the predicted and experimental KIE is not ideal. So, they performed molecular dynamics computations with the ONIOM approach using M06-2x/6-31G* for Reaction 1 and 53 THF molecules treated at PM3. 360 trajectories were begun in the region of the first transition state (TS1), and they can be organized into 4 groups. The first group (128 trajectories) are reactions that produce product. The second group (76 cases) form the C-C bond but then it ruptures and returns to reactant. The third group (82 cases) have an immediate recrossing back to reactant, and the last group (16 cases) takes product back to the first TS and then returns to product. The predicted KIE using this weighted MD results gives values in outstanding agreement with the experiments.

Of the first group, about 50% pass from TS1 to TS2 in less than 150 fs, or in other words look like a concerted path. But a good number of trajectories reside in the betaine region for 1-2 ps.

In contrast, trajectories initiated from the betaine with equilibrated THF molecules indicate a median of 600 ps to travel from TS1 to TS2 and do not resemble a concerted path.

They argue that this bimodal distribution is in part associated with a solvent effect. When the first TS is crossed the solvent molecules are not equilibrated about the solute, and 10-20% of the trajectories immediately pass through the betaine region due to “dynamic matching” where the entering motion matches with exiting over the second transition state. The longer trajectories result from improper dynamic matching, but faster motion in the solute than motion amongst the solvent needed to stabilize the betaine. So, not only do we need to be concerned about dynamic effects involving the reactants, we need to be concerned about dynamics associated with the solvent too!

References

(1) Chen, Z.; Nieves-Quinones, Y.; Waas, J. R.; Singleton, D. A. "Isotope Effects, Dynamic Matching, and Solvent Dynamics in a Wittig Reaction. Betaines as Bypassed Intermediates," J. Am. Chem. Soc. 2014, 136, 13122-13125, DOI: 10.1021/ja506497b.

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

Energy-conserving ab initio molecular dynamics

Contributed by David Bowler
Reposted from Atomistic Computer Simulations with permission

Energy conservation in molecular dynamics is an important test of the propagators and parameters being used.  In first principles Born-Oppenheimer molecular dynamics (BOMD) this is complicated because of the need for self-consistent iteration to a consistent ground state charge density and potential.  The self-consistency cycle introduces errors since it is not possible to achieve a perfect ground state.  Moreover, to improve efficiency, the starting charge density is often predicted from previous charge densities; while this can results in impressive speed gains, it breaks the time-reversibility of the dynamics and will introduce an unavoidable energy drift.

This blog will discuss a set of papers[1-4] from Anders Niklasson in Los Alamos, which seek to remove this problem, and which have also shown a way to avoid the self-consistency cycle entirely.  The work is perhaps less well-known that it should be, and offers an excellent route to efficiency.  The recent linear scaling DFT molecular dynamics that we have published (going to over 32,000 atoms with DFT)[5] used the time-reversible formulation to conserve energy.

The time reversibility is built on an extended Lagrangian (and is known as XL-BOMD) which introduces a new dynamical variable, or set of variables, for the electronic degrees of freedom (whether charge density or density matrix) alongside the nuclear degrees of freedom, in a trick reminiscent of Car-Parrinello MD[6].  This degree of freedom, written below as P, evolves harmonically in time, depending on a frequency (set by the user) and the difference to the ground state charge density (or whatever is being used).

This degree of freedom, P, can be propagated in a time-reversible manner with an accuracy that only depends on the timestep.  The BOMD can thus be made time-reversible by starting the self-consistency cycle from P – resulting in stability and energy conservation in MD over long timescales.  On its own, this is a very valuable advance.

However, in a set of papers[2-4], Niklasson takes the idea further.  He suggests that, given the right functional for the energy, the self-consistency cycle can be removed.  This involves performing a single diagonalisation/minimisation without self-consistent update, and using P as the input charge density.  The scheme has been developed, with the latest paper[4] generalising the earlier work, and showing that the approach is applicable to systems regardless of gap, provided that an adiabatic separation between electronic and nuclear degrees of freedom can be made.

This work provides techniques to accelerate ab initio MD significantly; we have found that, in combination with linear scaling DFT, we can access picoseconds of MD for tens of thousands of atoms even with a simple initial implementation.  The two approaches together promise to make large-scale, long time accurate MD significantly more affordable.

[1] Phys. Rev. Lett. 100, 123004 (2008) DOI:10.1103/PhysRevLett.100.123004

[2] J. Chem. Phys. 139, 214102 (2013) DOI:10.1063/1.4834015

[3] J. Chem. Phys. 140, 044117 (2014) DOI:10.1063/1.4862907

[4] J. Chem. Phys. 140, 164123 (2014) DOI:10.1063/1.4898803

[5] arxiv 1409.6085 and J. Chem. Theor. Comput. in press (2014) DOI:10.1021/ct500847y

Solution-Phase Dimerization of an Oblong Shape-Persistent Macrocycle

Chu, M.; Scioneaux, A. N.; Hartley, C. S.  J. Org. Chem. 2014, 79, 9009–9017
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

π-π-stacking has been a major theme of my blog, and is discussed in Chapter 3.5.4 in the Second Edition of my book. Most examples involved molecules that are nearly circular (like benzene or triphenylene). Hartley and co-workers discuss the π-π-stacking of the oblong molecule 1, comparing its experimental features with computed features of the model compound 2.¹

The key spectroscopic feature associated with assembly of 1 are the changes in the ¹H chemical shifts with increasing concentration. For example, the chemical shift of the three protons on the triphenylene unit shift upfield by 0.30 to 0.66 ppm as the concentration increases from 10^-5 to 10^-2 M.

To see if these NMR shift changes are due to association of 1, they employed a computational approach. First they optimized the structure of model compound 2 at B3LYP/6-31G(d) (Shown in Figure 1a). Then using this fixed geometry, they computed the ¹H chemical shifts of the dimer of 2. They explored the stacking distance (ranging from 3.2 to 4.0 Å along with varying the displacement of the two molecules along the major axis from 0.0 to 6.0 Å, finding the best fit to the chemical shifts with a separation of 3.6 Å and a displacement along the major axis of 3.5 Å. Using these two fixed values, they explored displacement of the molecules along the minor axis, along with rotation of the two molecules. The best fit to the experimental chemical shifts was with a displacement of 0.5 Å along the short axis and no rotation. This structure is shown in Figure 1b, with a RMS error of only 0.09 ppm from experiment. Models of the trimer show poorer fit to the experimental data.

(a) 2

(b) 2 dimer

Figure 1. B3LYP/6-31G(d) (a) optimized structure of 2 and the (b) structure of the best fit of the dimer of2. (As always, clicking on these images will allow you to manipulate the 3-D structure using JMol – highly recommended for the dimer.)

Using some smaller models and the B97-D functional, they argue that the displacement, which is substantially larger than the displacement found in stacked triphenylene, results from the need to minimize the steric interactions between the alkoxyl chains.

References

(1) Chu, M.; Scioneaux, A. N.; Hartley, C. S. "Solution-Phase Dimerization of an Oblong Shape-Persistent Macrocycle," J. Org. Chem. 2014, 79, 9009–9017; DOI: 10.1021/jo501260c.

InChIs:

1: InChI=1S/C116H148O10/c1-11-21-31-41-59-117-95-71-87-51-55-91-75-97-99(103-81-111(121-63-45-35-25-15-5)115(125-67-49-39-29-19-9)85-107(103)105-83-113(123-65-47-37-27-17-7)109(79-101(97)105)119-61-43-33-23-13-3)77-93(91)57-53-89-70-90(74-96(73-89)118-60-42-32-22-12-2)54-58-94-78-100-98(76-92(94)56-52-88(69-87)72-95)102-80-110(120-62-44-34-24-14-4)114(124-66-48-38-28-18-8)84-106(102)108-86-116(126-68-50-40-30-20-10)112(82-104(100)108)122-64-46-36-26-16-6/h69-86H,11-50,59-68H2,1-10H3
InChIKey=UUKOQEVXPHHWBJ-UHFFFAOYSA-N

2: InChI=1S/C76H68O10/c1-11-77-55-31-47-21-25-51-35-57-59(63-41-71(81-15-5)75(85-19-9)45-67(63)65-43-73(83-17-7)69(79-13-3)39-61(57)65)37-53(51)27-23-49-30-50(34-56(33-49)78-12-2)24-28-54-38-60-58(36-52(54)26-22-48(29-47)32-55)62-40-70(80-14-4)74(84-18-8)44-66(62)68-46-76(86-20-10)72(82-16-6)42-64(60)68/h29-46H,11-20H2,1-10H3
InChIKey=NGTKTKCZPONPTB-UHFFFAOYSA-N

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

Theoretical Foundation for the Presence of Oxacarbenium Ions in Chemical Glycoside Synthesis

Hosoya, T.; Takano, T.; Kosma, P.; Rosenau, T.  J. Org. Chem. 2014, 79, 7889-7894

Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

Comparing S_N2 and S_N1 reactions using computational methods is often quite difficult. The problem is that the heterolytic cleavage in the S_N1 reaction leads to an ion pair, and in the gas phase this is a highly endothermic process. Optimization of the ion pair in the gas phase invariably leads to recombination. This is disturbingly the result even when one uses PCM to mimic the solvent, which one might have hoped would be sufficient to stabilize the ions.

The computational study of the glycoside cleavage by Hosoya and colleagues offers some guidance towards dealing with this problem.¹ They examined the cleavage of triflate from 2,3,4,6-tetra-O-methyl-α-D-glucopyranosyl triflate 1.

Benchmarking the dissociation energy for the cleavage of 1 and considering computational performance, they settled on M06-2X/BS-III//M06-2X/BS-I, where BS-III is aug-cc-pVTZ basis set for O, F, and Cl and cc-pVTZ for H, C, and S and BS-I is 6-31G(d,p) basis sets were employed for H, C, and S, and 6-31+G(d) for O, F, and Cl. Solvent (dichloromethane) was included through PCM.

Attempted optimization of the contact ion pair formed from cleavage of 1 invariably led back to the covalent bound 1. PCM is not capable of properly stabilizing these types of ions in proximity. To solve this problem, they incorporated four explicit dichloromethane molecules. A minor drawback to their approach is that they did not sample much of configuration space and so their best geometries may not be the lowest energy configurations. Nonetheless, with four solvent molecules, they were able to locate contact ion pairs and solvent-separated ion pairs. Representative examples are shown in Figure 1. This method of explicit incorporation of a few solvent molecules seems to be the direction we must take to treat ions or even highly polar molecules in solution.

1 0.0	1-CIP 8.5
1-SSIPa 8.4	1-SSIPb 11.1

Figure 1. Representative examples of microsolvated 1, its contact ion pair (CIP) and solvent separated ion pair (SSIP) computed at M06-2X/BS-III//M06-2X/BS-I, and relative energies (kcal mol^-1)

References

(1) Hosoya, T.; Takano, T.; Kosma, P.; Rosenau, T. "Theoretical Foundation for the Presence of Oxacarbenium Ions in Chemical Glycoside Synthesis," J. Org. Chem. 2014, 79, 7889-7894, DOI:10.1021/jo501012s.

InChIs

1: InChI=1S/C11H19F3O8S/c1-17-5-6-7(18-2)8(19-3)9(20-4)10(21-6)22-23(15,16)11(12,13)14/h6-10H,5H2,1-4H3/t6-,7-,8+,9-,10-/m1/s1
InChIKey=RPZNYYCDDYUPJR-HOTMZDKISA-N

Update: be sure to read the interesting discussion associated with this post on Computational Organic Chemistry.


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