Potential Energy Surface-Based Automatic Deduction of

Article pubs.acs.org/JCTC

Potential Energy Surface-Based Automatic Deduction of Conformational Transition Networks and Its Application on Quantum Mechanical Landscapes of D‑Glucose Conformers Hiroko Satoh,*,†,‡,§,∥ Tomohiro Oda,⊥ Kumiyo Nakakoji,# Takeaki Uno,§ Hiroaki Tanaka,○ Satoru Iwata,○ and Koichi Ohno∥,∇ †

Research Organization of Information and Systems (ROIS), Tokyo 105-0001, Japan Department of Chemistry, University of Zurich, 8057 Zurich, Switzerland § National Institute of Informatics (NII), Tokyo 101-8430, Japan ∥ Institute for Quantum Chemical Exploration (IQCE), Tokyo 108-0022, Japan ⊥ Software Research Associates Inc., Tokyo 171−8513, Japan # Center for the Promotion of Interdisciplinary Education and Research, Kyoto University, Kyoto 606−8501, Japan ○ Department of Mathematical Informatics, University of Tokyo, Tokyo 113−8654, Japan ∇ Department of Chemistry, Graduate School of Science, Tohoku University, Sendai 980−8578, Japan ‡

S Supporting Information *

ABSTRACT: This paper describes our approach that is built upon the potential energy surface (PES)-based conformational analysis. This approach automatically deduces a conformational transition network, called a conformational reaction route map (r-map), by using the Scaled Hypersphere Search of the Anharmonic Downward Distortion Following method (SHSADDF). The PES-based conformational search has been achieved by using large ADDF, which makes it possible to trace only low transition state (TS) barriers while restraining bond lengths and structures with high free energy. It automatically performs sampling the minima and TS structures by simply taking into account the mathematical feature of PES without requiring any a priori specification of variable internal coordinates. An obtained r-map is composed of equilibrium (EQ) conformers connected by reaction routes via TS conformers, where all of the reaction routes are already confirmed during the process of the deduction using the intrinsic reaction coordinate (IRC) method. The postcalculation analysis of the deduced rmap is interactively carried out using the RMapViewer software we have developed. This paper presents computational details of the PES-based conformational analysis and its application to D-glucose. The calculations have been performed for an isolated glucose molecule in the gas phase at the RHF/6-31G level. The obtained conformational r-map for α-D-glucose is composed of 201 EQ and 435 TS conformers and that for β-D-glucose is composed of 202 EQ and 371 TS conformers. For the postcalculation analysis of the conformational r-maps by using the RMapViewer software program we have found multiple minimum energy paths (MEPs) between global minima of 1C4 and 4C1 chair conformations. The analysis using RMapViewer allows us to confirm the thermodynamic and kinetic predominance of 4C1 conformations; that is, the potential energy of the global minimum of 4C1 is lower than that of 1C4 (thermodynamic predominance) and that the highest energy of those of all the TS structures along a route from 4C1 to 1C4 is lower than that of 1C4 to 4C1 (kinetic predominance).

1. INTRODUCTION

however, has been found very difficult or sometimes impossible even with those advanced experimental techniques. Conformational sampling and the global minimum problem have also been a major topic in computational chemistry.1 The goal of conformational analysis is to fully cover the global conformational space and determine all of the minima and

Conformational change of chemical structures is an important factor in determining their chemical properties and behaviors of molecular systems especially for those with high flexibility. Existing experimental methods to access three-dimensional characteristics of molecules include NMR spectroscopy (e.g., Nuclear Overhauser Effect (NOE), J-values, and chemical shifts), X-ray crystallography, vibrational spectroscopy, and circular dichroism spectroscopy. Capturing all conformers, © XXXX American Chemical Society

Received: April 29, 2016

A

DOI: 10.1021/acs.jctc.6b00439 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

information. These methods are more practical for a local conformational search. The random search and the MD simulation techniques can locate the minima, but neither of these techniques provides information about the pathways between the minima or among their transitions. Therefore, any of those methods requires postcalculation or postsimulation data analyses to locate the pathways between the stationary points on the obtained PES or FES. A nudged elastic band (NEB) method and its variations8−10 are commonly used for this purpose. Kirkov and Karplus applied clustering and statistical methods to reveal the peptide-folding map on FES obtained by MD simulations.11 Noé and Fischer reported an approach based on discrete mathematical methods to construct transition networks of macromolecules from MD trajectories.12 Carr and Wales introduced a method based on a graph theoretical algorithm with a doubly NEB method13−15 for finding pathways to deduce energy landscapes from a collection of the minima. Grouleff and Jensen performed a conformational search of peptides16 by locating the minima with MC calculations followed by guessing TS conformers on the force field energy surface by using three different methods: quadratic synchronous transit (QST),17,18 one type of the NEB method, and Scaled Hypersphere Search (SHS).19−22 Many of those postsimulation or postcalculation data analyses based on PES or FES, such as the series of NEB methods, can be classified into a PES- or FES-based data analysis. In contrast, there is a method that can be classified into a PES-based search. This is the low-mode search (LMOD) method,23 which generates structures based on PES by searching the low-frequency modes systematically or stochastically. The LMOD method makes it possible to efficiently locate minima and find even all the minima for a small molecular system. The LMOD method also needs the postcalculation data analyses for the obtained minima to locate transitions and pathways between them. Our approach we introduce in this article is built upon the PES-based conformational search and data analysis. We use the SHS of the Anharmonic Downward Distortion Following method (SHS-ADDF)19−22 to explore a global conformational space while focusing simply on the mathematical features of the potential energy hypersurface. The exploration almost fully covers the global conformational space. The outcome is a conformational transition network, which we call a conformational reaction route map (r-map), composed of nearly all the equilibrium (EQ) and TS conformers as well as pathways confirmed with the intrinsic reaction coordinate (IRC) method. Our method locates all of the EQ and TS structures as well as the pathways during the process of the deduction of an r-map. All these processes are automatically performed without requiring any a priori specification of variable internal coordinates. It is also notable that the result of the exploration is independent of the initial conformer. The exploration time has been significantly reduced from that required for a systematic search by searching rational structures along a PES. Our postcalculation data analysis is carried out for the obtained r-map by using the RMapViewer software we have developed to interactively visualize and analyze energy and structure data as well as the pathways in the transition network. Functions implemented in RMapViewer include to visualize a transition network or an energy landscape and to search the minimum energy paths (MEPs) based on the energy values of the EQ and TS structures with the pathway information. We

transition states (TS) as well as the pathways connecting them. It would be ideal if the calculations are performed for a molecular system in such a condition that is similar to nature (e.g., in solvent) using a high level theory with low computational cost. Computational methods to deal with this issue have been found powerful and useful and often been employed complimentarily with experimental measurements,2 and yet exploring the global conformational space to reveal global conformational networks containing all the minima and their transitions with minimum energy pathways is still a highly challenging problem. This paper describes our approach that makes it possible to nearly fully explore the global conformational space with low computational cost. Existing computational techniques can be classified based on how to perform the sampling and which type of theory is used. The approach we present here is a type of potential energy surface (PES)-based analysis based on the quantum mechanical (QM) theory. The primary difference from the existing methods includes the data contents obtained in the course of the exploration and the role of the postcalculation data analysis. The full coverage of the global conformational space can be theoretically achieved by a systematic sampling for all of the changeable internal coordinates. This geometry-basis comprehensive sampling, however, often quickly encounters a combinatorial explosion problem.2,3 Thus, a global systematic geometry search is sometimes regarded as theoretically possible but practically impossible for a large number of molecules and as feasible only for very small or rigid systems with little geometrical freedom. Geometry-based systematic sampling has thereby been primarily employed for local conformational search by focusing on selected bonds. The geometry-basis random search4,5 is often used for sampling the minima from a large conformational space with low computational cost. The Monte Carlo (MC) and the stochastic methods for the random search work well for efficient sampling of various minimum conformers. Those methods are sometimes used with a molecular dynamics (MD) technique. Molecular dynamics (MD) simulation is effective to explore conformational diversity. The MD techniques allow us to access thermodynamic behaviors of various types of molecular systems in conditions similar to nature. An MD simulation results in a global conformational search if an equilibrium simulation is fully completed. Such a full simulation, however, usually requires immensely long simulation time to complete all of the possible events. There is another problem in fully covering the global conformational space by MD simulations as pointed out in the textbook by Jensen.2 That is, the primary disadvantage of MD is the inability to overcome barriers larger than the internal energy determined by the simulation temperature. Sampling techniques, such as local elevation6 and metadynamics,7 have been developed to deal with this problem. These techniques are known to be effective to accelerate the simulation of a rare event and to reconstruct a free-energy surface (FES) for the event. The FES landscape reconstructed via local elevation or metadynamics methods contains information about all the minima as well as their transitions in the focused event. The practical number of collective variables is, however, up to three or four because a simulation with a larger number of collective variables demands exponentially large computational resource. The analysis of the results becomes extremely challenging to obtain any useful B


Article


another EQ structures. These functions are also used for analyzing conformational r-maps.

have applied this approach to a conformational analysis of Dglucose. In what follows, we first describe the detailed method of the PES-based conformational search using SHS-ADDF and the postcalculation data analysis using RMapViewer. We then present its application to deducing conformational r-maps of Dglucose. We describe the analysis of the conformational r-maps by demonstrating the use of RMapViewer with a particular focus on our findings of the minimum energy paths (MEPs) of the flipping reactions between 4C1 and 1C4 ring conformations. Some of our colleagues have recently published a preliminary letter briefly reporting an application of the SHS to a conformational search of L-serine.24 This article provides more details about the computational methods of the PESbased conformational search, as well as demonstrates an extensive conformational analysis using RMapViewer, which helps researchers to effectively derive valuable information about D-glucose from the obtained conformational r-maps.

3. APPLICATION TARGET We have chosen a simple application target, D-glucose, in this study to demonstrate the efficiency of our conformational analysis based on potential energy hypersurface. The reason we have chosen a monosaccharide is that the conformational behavior of a monosaccharide has been a primary interest in carbohydrate chemistry. It continues to be intensively studied today ever since the early study on conformational free energy of carbohydrate was presented.37 The conformational analysis of monosaccharides provides fundamental data to investigate reactivity and functions of carbohydrates in both chemical and biological systems.38 For instance, the relationships between sugar conformations and their reactivity and selectivity in glycosylation reactions have been intensively studied in carbohydrate chemistry.39−42 We have conducted the exploration in the QM level because they must be treated at the QM level to take into account the stereoelectronic effects, such as an anomeric effect.43,44 The stereoelectronic effect theory is often used to discuss the stabilization and destabilization of conformers of saccharides. This specific feature of saccharides poses a challenge for the conformational analysis of this class of molecules difficult. Glucose is a simple and quintessential monosaccharide. QM conformational studies on D-glucose have been intensively performed since the late 1970s.45−55 QM studies toward energy landscapes of the global conformation space have been reported at various levels of theories especially over the past decade.48,49,51,53−55 Ionescu et al. revealed conformational inversions between 4C1 and 1C4 conformations of D-glucose by using static and dynamical density functional calculations.48 The first trial of deducing a QM-based energy landscape of monosaccharide conformers was conducted by using the first principle Car−Parrinello molecular dynamics simulations.49 Construction of QM-based kinetic conformational landscapes for five kinds of monosaccharide has recently been completed by comprehensive generations of geometries of an isolated molecule in the gas phase using MM computational tools, followed by QM geometry optimizations to locate EQ and TS structures, and allocations of the reaction routes between them with IRC.55 We have deduced a conformational landscape for an isolated molecule in the gas phase. We have chosen such an ideal model as the first application, because it is still essentially important to provide ample basic information about geometric and electronic structures. It would stimulate ideas about new synthetic and molecular design and for solving the mechanisms of biochemical systems.

2. METHODS 2.1. SHS-ADDF and GRRM. SHS-ADDF, which serves as the computational and theoretical foundation for the PES-based conformational search, was originally developed by Ohno and Maeda.19−22 SHS-ADDF was initially designed to explore all possible chemical reaction events associated with an input chemical formula. SHS-ADDF has been implemented as a part of the GRRM program package. The output of SHS-ADDF is a global reaction route map (r-map) consisting of all of the chemical reaction events composed of EQ structures and dissociation channels (DCs) connected by reaction routes via TS structures. All of the reaction routes are confirmed using the IRC method. Existing studies have applied SHS to investigate mechanistic details of various types of chemical reactions. Ohno and Maeda discussed the applications of SHS in its early stages (2004−2005).20−22,25−27 The latest reviews and in the references therein provide details of SHS and its applications.28,29 The core idea of the method we propose in this paper is the application of SHS-ADDF to conformational search by following only large ADD, by restricting bond lengths, and by cutting off structures with high free energy. This simple idea has made it possible to automatically perform sampling only rational structures in a conformational space resulting in lower cost of computational time and resources even at the quantum mechanical (QM) level. 2.2. RMapViewer. We have been studying the PES-based conformational analysis in the Maizo-chemistry project30 for the last several years. The Maizo-chemistry project aims at molecular and reaction discovery. We have applied SHS-ADDF to obtaining the global reaction route mappings by performing explorations without any restrictions31 and have discovered a new carbon family as Maizo-molecules.32−34 In the course of the Maizo-chemistry project, we have developed an application system, RMapViewer,35 for uses to interactively analyze r-maps. Since an r-map usually results in a large network structure,36 we believe that a visualization software tool to help researchers analyze the r-maps with such complex contents is essential. We have implemented several functions designed for the analysis of the r-maps on RMapViewer, including functions for searching for structures, for visualizing r-maps, molecular structures, and the molecular structure change along an energy profile, as well as for enumerating all of the possible routes from one EQ to

4. COMPUTATIONAL DETAILS Figure 1 illustrates the procedural steps of the PES-based conformational analysis using GRRM and RMapViewer. The study used SHS-ADDF implemented in GRRM1156 for deducing r-maps. The study then employed the RMapViewer 4.0 software for the postcalculation data analysis of the obtained r-maps. All of the QM calculations were performed using the restricted Hartree−Fock (RHF) method at the 6-31G level with the Gaussian 09 program package.57 We used a computer with Intel Xeon Processor E5-2667 v2 for the exploration with SHS-ADDF. C


Article


number without unit. The smaller the LADD value is chosen, the larger ADDF becomes, meaning to explore only reaction routes connected via lower TS barriers. Using small LADD values may result in missing conformers via high TS barriers, but using large LADD values may lead to deal with more bondbreaking channels. The goal here is to find a good tuning point that detects more conformers and less bond-breaking channels so that largely populated conformers can be effectively explored with relatively low cost of computational time and resources. Thus, the LADD option is effective to efficiently trace reaction routes via low TS barriers on a PES. Note that there is no guarantee that a conformational search can be carried out entirely if one uses the LADD option. Nevertheless, it can be guaranteed that largely populated conformers, which are mostly important, would all be efficiently obtained if one chooses an appropriate LADD value. This study set LADD = 5, which resulted in r-maps composed of mostly conformers and a few bond-breaking channels. The second GRRM option is called bond condition, which restrains bond lengths between atoms m and n, lbond(m,n), not to trace further if the molecule undergoes bond-breaking and/ or bond-formation. This study presumed that if the distance between atoms m and n in the initial geometry, lbond(m,n)initial, was less than or equal to 1.80 Å, a bond was assumed to exist between the pair. For each bonding pair, the maximum bond length, lbond(m,n)max was set as follows in Å

Figure 1. PES-based conformational analysis using GRRM and RMapViewer. All of the EQ and TS structures are conformers.

4.1. Preparation of the Initial Structures of Glucose. The initial structures of D-glucose (1 and 2 in Figure 2) were generated via full geometry optimization with the Opt = tight option in the gas phase for the anomeric isomers (anomers) α (1) and β (2). Optimizations were initiated from an ideal 4C1 chair conformation (1a and 2a). Note that the search result is independent of which initial structure is chosen; that is, the same r-map is derived starting from any initial conformer. This is so because SHS-ADDF traces all the barriers if any restrictions are specified. In this study SHS-ADDF traces all the barriers lower than the limitation, which is determined with the GRRM options not to break bonds (see Section 4.2). Therefore, all the conformers produced via bond rotations can be explored starting from any initial conformer.

lbond(m , n)max = lbond(m , n)initial × 1.20

During the calculation, when lbond(m,n) exceeds lbond(m,n)max, further exploration was not carried out beyond the geometry. The third GRRM option is NLowest. This option is used to restrain the SHS from searching the structures with high free energies. The NLowest is a number without unit. We set this value at 40 in this study. These options control the GRRM exploration not to trace beyond the restrictions. That is, if the size of ADD is smaller than the LADD restriction, the height of the TS barrier is estimated to be large based on the ADD model, and further trace toward the TS will not be performed. If an obtained EQ structure is with high free energy and/or if any of the bond length of it is found to be longer than the maximum bond length (lbond(m,n)max), further routes from the EQ structure will not be traced, either. Therefore, starting with an initial structure in high energy might cause a termination of the exploration before it is completed, even though the GRRM exploration is independent of its initial structure as mentioned above. In general, choosing an initial structure in low energy is safer. One more thing to note here is that if another conformer network exists beyond the restrictions, the network will not be explored. For instance, if a conformer produced via an anomerization reaction belongs to another isolated conformer network, neither of the conformer or the network will be explored. Such GRRM options together with the xyz coordinates of the initial structure are specified in the input file to GRRM. Full input files of the GRRM calculations are provided in the Supporting Information. For more detailed information about these GRRM options, see ref 58. 4.3. Refinement of the r-Maps. After the GRRM calculations, we examined if the obtained r-maps contain any reacted structures connected via bond-breaking channels by using the CAST (CAnonical representation of STereo-

Figure 2. Targets of conformational analysis. α-D-glucose (1) and its two types of chair conformations, 4C1 (1a) and 1C4 (1b). β-D-glucose (2) and its two types of chair conformations 4C1 (2a) and 1C4 (2b).

4.2. PES-Based Conformational Search Using GRRM. The procedure starts with PES-based conformational search with SHS-ADDF (Figure 1a). As described in Section 2, the PES-based conformational search is carried out by tracing only large ADD, restraining bond lengths, and cutting off structures with high free energy. The size of the ADDF can be tuned by using the large-ADDF (LADD) option of the GRRM program. The LADD value is a D


Article


Figure 3. A global view of the r-map for α-D-glucose (1) visualized with RMapViewer 4.0.

chemistry) line notation.59,60 Such reacted structures are then removed from the r-maps, which results in conformational rmaps, composed of only those EQ conformers that are connected by reaction routes via TS conformers. 4.4. Analysis of the r-Maps. We classified the EQ and TS conformers into the 38 types of conformers and located the global minima of a 4C1 conformation for the anomers of α (1amin) and β (2amin) and those of a 1C4 conformation for anomers α (1bmin) and β (2bmin). The classification conforms to the IUPAC nomenclature for pyranoses.61 In the IUPAC rule, five types of conformations are defined for six-membered rings: chair (C), boat (B), twist/skew-boat (S), half-chair (H), and envelope (E) conformations. Regarding pyranoses, the serial numbers of carbon atoms (Figure 2) are added to those conformer notations according to the relative positions compared to a reference plane, and 38 types of conformers are distinguished. These conformers are described numerically by using the Cremer-Pople puckering coordinates (CP coordinates), which is the polar coordinates system plotted on a sphere (CP sphere).62 We calculated the CP coordinates using the program developed by Whitfield et al.63 The last step of our conformational analysis is the postcalculation analysis of the conformational r-maps using RMapViewer 4.0 (Figure 1b). We retrieved all possible reaction routes from 1amin (4C1) to 1bmin (1C4), from 1bmin (1C4) to 1amin (4C1), from 2amin (4C1) to 2bmin (1C4), and from 2bmin (1C4) to 2amin (4C1). Once the starting node (i.e., a reactant) and the ending node (i.e., a product) are interactively specified

by clicking on the nodes displayed in the map, RMapViewer automatically enumerates all possible routes from the starting node to the ending node. The enumerated routes are sorted in the ascending order of the TSmax value, which is the highest of all the TS energies along one route. The sorting process also takes into account the number of hops, where one hop corresponds to a hill from one EQ to another EQ via one TS. As a result, the shortest routes with the smallest TSmax value are output to be the MEPs. One may also specify the maximum number of hops and routes to enumerate. We used default values: 8 hops and 100 routes in this study.

5. RESULTS AND DISCUSSIONS This section presents the results of applying PES-based conformational analysis to D-glucose and discusses the findings of the results. 5.1. Obtained Structures in r-Maps. The obtained conformational r-map for α-D-glucose (1) was composed of 201 EQ and 435 TS conformers. The r-map resulted from GRRM containing bond-breaking channels, which consist of reacted structures that underwent bond breaking/forming, was composed of 209 EQ, 443 TS, and 1 DC structures. The r-map was obtained via 287,686 force and 17,932 Hessian calculations. The consumed CPU time was 98.9 h. The bond-breaking channels (8 EQ, 8 TS, and 1 DC) include a C1−O5 endo-cyclic cleavage reaction, which is initiated by the H1 proton transfer to O5 to produce an aldohexose. Its TS energy was 222.9 kJ/mol. This relatively high TS energy must be due to the fact that any E


Article


Figure 4. EQ conformers (red circles) and TS conformers (blue squares) plotted in the Mercator projections of the Cremer-Pople puckering sphere (CP sphere).

for the α-anomer (1). The EQ conformers of the β-anomer (2) do the same but with one exception, adopting a 3 H 2 conformation. The TS conformers occupy a large conformational space with somewhat sparse regions. The EQ conformer highest in energy of the α-anomer (1) adopts a 1S3 conformation, 73.6 kJ/mol higher than that of the global minimum (1amin). The TS conformer lowest in energy adopts a 4C1 conformation, 12.6 kJ/mol higher than that of 1amin. The TS conformer highest in energy adopts a B3,O conformation, 94.6 kJ/mol higher than that of 1amin. For βanomer (2), the EQ conformer highest in energy adopts a 5S1 conformation, 72.2 kJ/mol higher than that of the global minimum (2amin). The TS conformer lowest in energy adopts a 4 C1 conformation, 15.1 kJ/mol higher than that of 2amin. The TS conformer highest in energy adopts an E1 conformation, 74.9 kJ/mol higher than that of 2amin. It is interesting that rather wide ranges of energy are found for geometries in the same conformational group for both EQ and TS conformers. For example, the energy range of the EQ conformers adopting 4C1 is about 50 kJ/mol for both anomers 1 and 2. This means that the rotation of the side-chains is an important factor in stabilizing and destabilizing the glucose structures. This may also infer that the formation of the sidechains plays an important role in the course of the

counterions, which stabilize the system, have not been considered in the calculations. Figure 3 shows a global view of the r-map of 1 visualized with RMapViewer. This is a twodimensional projection of a conformational transition network containing EQ and TS conformers connected with reaction pathways on hypersurface. The obtained conformational r-map for β-D-glucose (2) was composed of 202 EQ and 371 TS conformers. The r-map resulted from GRRM including reacted structures was composed of 205 EQ, 373 TS, and 1DC structures, which was obtained via 357,803 force and 14,418 Hessian calculations. The consumed CPU time was 95.2 h. The global minimum adopts a 4C1 conformation for both anomers α and β (1amin and 2amin, respectively). The 1amin with the axial hydroxyl group at the anomeric site is expected to be stable than 2amin because of the anomeric effect. The calculations showed this tendency that the 1amin is 10.1 kJ/ mol lower than 2amin in energy. In Figure 4, the distributions of conformers in the conformational r-maps are plotted in the Mercator projections of the CP sphere.62 In Figure 5, each conformer is plotted against the potential energy relative to 1amin or 2amin. In the landscapes reported by Beckham et al.,55 all the EQ conformers adopted chair (C), boat (B), or skew (S) conformations. Our calculations show the same tendency: all the EQs adopt chair (C), boat (B), or skew (S) conformations F


Article


Figure 5. EQ conformers (red bars) and TS conformers (blue stars) classified into 38 conformational types are plotted against potential energy.

with the pyranose conformational change of 4C1 → E3† → B3,O, as shown in Figure 8, which takes place after 3 hops where the rotation of the side-chains takes place at C1−O1, C2−O2, C3− O3, C4−O4, C5−C6, and C6−O6. The difference between groups-1 and -2 is the pyranose conformational change after the TSmax hop. The pyranose rings of group-1 undergo 2E → 5S1 → 5 H4 → 1C4 changes, while those of group-2 undergo 1S3 → 1,4B → OS2→ 3H2 → 1C4 changes. Two MEPs with 11 hops were found in the opposite direction: from 1bmin (1C4) to 1amin (4C1) (Figure 9). The trajectories and energy profiles are shown in Figures 9a and 9b, respectively. The TSmax hill is located at the first hop for both MEPs, with 15.6 kJ/mol, where the rotation of C5−C6 and C3− O3 simultaneously takes place as the pyranose ring keeps the form as a 1C4 conformation, as shown in Figure 10. The torsion angles γ1 and γ2, corresponding to the rotation at C5−C6 and C3−O3, respectively, are also listed in Figure 10. These angles show that both of the conformations associated with the C5−C6 and C3−O3 bonds at TSmax adopt an eclipsed conformation, whereas those in the neighboring EQ structures adopt a more stable staggered conformation. After the TSmax hill, via some hops with the rotations of the side-chains, the pyranose ring undergoes 3H2 → OS2 → B2,5→ 1,4B → 1S3→ 4H3 → 4C1 changes. For β-anomer (2) from 2amin (4C1) to 2bmin (1C4), three MEPs were found with 8 hops (Figure 11). The trajectories and energy profiles are shown in Figures 11a and 11b, respectively.

conformational change of a pyranose ring or when undergoing chemical reactions. 5.2. Finding MEPs between 4C1 and 1C4 conformers. The MEP search on the r-maps was carried out using RMapViewer 4.0, which has been designed to help researchers in exploring the space of the r-maps. A snapshot of the analysis of α-anomer (1) is shown in Figure 6. This is a step to analyze the second route out of the enumerated routes from 1amin (4C1) (a square in the single blue line in the map displayed in the top left frame) to 1bmin (1C4) (a square in the double blue lines in the map). The second route is highlighted in the map, and the list of the enumerated routes is displayed in the bottom frame. The reaction movie (shown in the top center of Figure 6) associated with the reaction profile plot (in the top right) helped us visually analyze the conformational change along the selected route. For α-anomer (1) from 1amin (4C1) to 1bmin (1C4), 10 MEPs with 7 hops were found (Figure 7). Regarding the similarity of the conformational change of the pyranose ring, they are classified into two groups: paths 1−5 (group-1) and paths 6− 10 (group-2). The highlighted route in Figure 6 corresponds to path 4. The trajectories and energy profiles of group-1 are shown in Figures 7a and 7b, respectively. Those of group-2 are shown in Figures 7c and 7d. The trajectories are plotted in the Mercator projection of the CP sphere. The TSmax hills are observed at the fourth hop for all 10 paths, marked with a large blue circle in Figure 7. The height of the hill is 33.7 kJ/mol, G


Article


Figure 6. RMapViewer 4.0 in the analysis of MEPs from 1amin (4C1) to 1bmin (1C4). Top left: a map view. Bottom frame: a list of the enumerated routes from the starting node (a square in blue single line in the map view) to the ending node (a square in double blue lines in the map view). Top center: molecular view and movie. Top right: reaction profile plot. The second route out of the enumerated routes is selected to analyze, which is highlighted in the map and the route list. The conformational and energy change is visually analyzed using the movie associated with the reaction profile.

its derivatives. For example, this information would be useful to modify structures of glycosyl donors for high selectivity and/or reactivity by employing appropriate protecting groups at sidechains that can control the conformational transformation of pyranose ring.39 The influence of the side-chain conformers on the pyranose ring conformers is interesting to study, although we would need more studies to derive a definite conclusion in addition to the currently obtained results. Glucose is assumed to favorably adopt a 4C1 conformation rather than a 1C4 conformation to avoid the steric repulsion.38 The calculations showed this tendency in the comparison between 4C1 and 1C4 conformations (Figure 5). Regarding the global minima of 4C1 or 1C4 conformations, 1amin (4C1) is 19.1 kJ/mol stable than 1bmin (1C4), and 2amin (4C1) is 18.0 kJ/mol stable than 2b min ( 1 C 4 ) in energy. Furthermore, the predominance of 4C1 was supported by the MEP analysis: For the α-anomer (1), the TSmax barrier from 4C1 to 1C4 is 18.1 kJ/mol higher than the TSmax barrier in the opposite direction. For the β-anomer (2), the TSmax barrier from 4C1 to 1C4 is 3.7 kJ/mol higher than the TSmax barrier in the opposite direction. Although the predominance of sugar conformations is usually discussed only with the stability of equilibrium structures, meaning, only with the thermodynamic factor, however, our calculations revealed that adopting 4C1 conformations is not only thermodynamically but also kinetically favorable in glucose as discussed by Ionescu et al.48 Furthermore, while the calculations show the thermodynamic predominance of 4C1 for the α-anomer higher than the β-anomer as mentioned above, the MEP analysis shows also the kinetic predominance of 4C1 of the α-anomer higher than the β-anomer.

The TSmax hill is located at the second hop for all three MEPs, with 16.2 kJ/mol, where the pyranose ring changed the conformation to 4C1 → 2H3† → B3,O, as shown in Figure 12. The change in the pyranose ring at the TSmax hill is very similar to the case of α-anomer (1). After the TSmax hill, there are two branches. One is B3,O → 1S3 → 1,4B → 1H2 → 1C4, and the other is B3,O → 1S3 → 1,4B → 1S5 → OS2 → 3H2 → 1C4. From 2bmin (1C4) to 2amin (4C1), 4 MEPs with 15 hops were found (Figure 13). The trajectories and energy profiles are shown in Figures 13a and 13b, respectively. All of the paths show the same change of the pyranose ring conformation. The TSmax hill is located at the second hop for all four MEPs, with 12.5 kJ/mol, where only the rotation of C5−C6 takes place as the pyranose ring keeps the form of an 1C4 conformation, as shown in Figure 14. The torsion angles γ1 corresponding to the rotation at C5−C6 are also listed in Figure 14. These angles show that the conformation associated with the C5−C6 bond at TSmax adopts an eclipsed conformation, whereas those in the neighboring EQ structures adopt a more stable staggered conformation. At the first hop, another side chain rotation at C6−O6 takes place. After the TSmax hill, via some hops with the rotations of the side-chains, the pyranose ring undergoes 3H2 → O S2 → B2,5→ 1S5→ B3,O→ 2H3 → 4C1 changes. Those MEP analyses have indicated that not only the ring conformational changes but also the side-chain rotations play an important role in stabilizing or destabilizing the molecular system. The molecules seem to prepare for the transformation of ring conformers by rotating the side-chains, or the rotations of the side-chains seem to induce the ring conformational change. These behaviors are informative for us in the analysis of the synthetic or biochemical transformations of glucose and/or H


Article


Figure 7. MEPs from 4C1 (1amin) to 1C4 (1bmin) for α-D-glucose (1). (a) Plots in Mercator projections of CP sphere for paths 1−5. (b) Energy profiles of paths 1−5. (c) Plot in Mercator projections of CP sphere for paths 6−10. (d) Energy profiles of paths 6−10. I


Article


Figure 8. Geometries of 1amin, TSmax and neighboring EQ conformers of MEPs from 1amin (4C1) to 1bmin (1C4).

Figure 9. MEPs from 1C4 (1bmin) to 4C1 (1amin) for α-D-glucose (1). (a) Plots in Mercator projections of CP sphere. (b) Energy profiles.

time was 236.0 h. Note that many more reacted structures were explored than that with LADD = 5. Newly explored conformers in this calculation were OE, 2H1, E2, and 3H4 conformations, which are all the TS conformers. Some new TS and EQ conformers belonging to the conformational types explored with LADD = 5 were also found mostly in high energy levels, especially for E 5, E 1, 4 H 3, 1,4B, 1 S5, B3,O , 1 S3, and 5H 4 conformations. Most of these conformers were shown as high TS conformers in the landscape for α-D-glucose reported by Beckham et al.55 Although the QM level in this study is different from their calculations, the tendency of our results seems to be conforming to their results. The obtained conformer types plotted against the potential energy relative to 1amin are provided in the Supporting Information.

Figure 10. Geometries of 1bmin, TSmax, and its neighboring EQ conformers of MEPs from 1bmin (1C4) to 1amin (4C1). At the hop of TSmax, no ring conformational change but bond rotations at C5−C6 and C3−O3 simultaneously take place. The torsion angles corresponding to these rotations are listed.

6. INFLUENCE OF THE QM LEVEL The aforementioned discussions are based on the calculations conducted at the RHF/6-31G level. We used this moderate level in our study because the primary purpose of this article is to demonstrate the effectiveness of the proposed method. Further calculations at a higher level may be necessary for more detailed investigations on specific problems. In what follows, we discuss the influence of the QM levels on the PES-based conformational search using SHS-ADDF. We chose RHF/6-31G as a moderate level based on the previous discussions by Maeda and Ohno on the influence of the QM level on the exploration with SHS-ADDF in calculating H2O octamers.58 They showed that the geometries in lower potential energy are mostly the same between RHF/6-31G and B3LYP/6-311+G(d,p) levels, although the number of EQ

5.3. Conformational Search with a Larger LADD Value. To test the influence of the LADD value, we have conducted conformational search for α-D-glucose (1) with the same parameters and the initial structure but with a larger LADD = 10, which widens the PES space to explore for higher TS barriers than the case discussed above with LADD = 5. As expected, we obtained a larger conformational r-map composed of 222 EQ and 644 TS conformers. The r-map resulted from GRRM including reacted structures was composed of 400 EQ, 763 TS, and 5 DC structures, which was obtained via 1,102,477 force and 56,466 Hessian calculations. The consumed CPU J


Article


Figure 11. MEPs from 4C1 (2amin) to 1C4 (2bmin) for β-D-glucose (2). (a) Plots in Mercator projections of CP sphere. (b) Energy profiles.

Figure 12. Geometries of 2amin, TSmax, and its neighboring EQ conformers of MEPs from 2amin (4C1) to 2bmin (1C4).

Figure 13. MEPs from 1C4 (2bmin) to 4C1 (2amin) for β-D-glucose (2). (a) Plots in Mercator projections of CP sphere. (b) Energy profiles. K


Article


Figure 14. Geometries of 2bmin, TSmax, and its neighboring EQ conformers of MEPs from 2bmin (1C4) to 2amin (4C1). At the hop of TSmax, no ring conformational change but bond rotation at C5−C6 takes place. The torsion angles corresponding to this rotation are listed.

conformational search by comparing with the computationally obtained conformers in low energy giving less than 10 kJ/mol relative to 1amin for the α-anomer (1) and 2amin for the βanomer (2) at the RHF/6-31G level. The PES-based conformational search for 1 found five conformers (1amin and 1a2−1a5) to be in low energy giving less than 10 kJ/mol relative to 1amin. Figure 16 shows these geometries with the intramolecular O···H distances shorter than 3.0 Å. The computational search results include all of the experimentally observed four conformers labeled as G−g+/cc/ t, G+g−/cc/t, Tg+/cc/t, and G−g+/cl/g− corresponding to 1a3, 1a4, 1a2, and 1a5, respectively. The O···H distances are in good agreement with those reported by Alonso et al.64 However, the computational ordering in energy (1amin → 1a2 → 1a3 → 1a4 → 1a5) is not in good agreement with the experimental observations (1a3 and 1a4 → 1a2 and 1a5), and the conformation of the computational minimum 1amin is not included in the experimental report. We also performed single point energy calculations at the B3LYP/6-31G and MP2/6-311++G(d,p) levels (Table 1 footnotes c and e, respectively) for the RHF/6-31G geometries to investigate the influence of the QM level. In Figure 16, the obtained energies relative to 1amin at RHF/6-31G, B3LYP/631G, and MP2/6-311++G(d,p) are given from the left side in the parentheses under each of the geometries. The ordering at the B3LYP/6-31G level (1a2 → 1amin → 1a3 → 1a4 → 1a5) still is not in good agreement with the experiments, but excellent agreement was found at the MP2/6-311++G(d,p) level (1a3 → 1a4 → 1a2 → 1a5 → 1amin). We performed also geometry optimization calculations at these levels (Table 1 footnotes d and f, respectively) and found only very slight changes in energies and geometries that hold the same conformation types even after the geometry optimizations. For the β-anomer (2), three conformers (2amin and 2a2−2a4) were found to be in low energy giving less than 10 kJ/mol relative to 2amin by the PES-based conformational search. Figure 17 shows these geometries with the intramolecular O··· H distances shorter than 3.0 Å. These three conformer types are exactly the same as those of the three experimentally observed conformers labeled as G−g+/cc/t, G+g−/cc/t, and Tg+/cc/t corresponding to 2amin, 2a3, and 2a2, respectively. The O···H distances are in good agreement with those reported by Alonso et al.64 Only one difference between computational and experimental results is the ordering of these conformers in energy at 2a2 and 2a3. Then, we performed at the B3LYP/6-31G and MP2/6-311+ +G(d,p) levels single point energy calculations for the RHF/631G geometries and geometry optimization calculations for these three conformers. The single point energies relative to 2amin at the B3LYP/6-31G and MP2/6-311++G(d,p) levels are shown in Table 1 footnotes c and e, respectively. These energies are given under the geometries in Figure 17 in the

structures was reduced from 181 (RHF/6-31G) to 164 (B3LYP/6-311+G(d,p). We assumed that the same trend would be observed in this study of the application to conformational search. We have verified this assumption by comparing the results at the RHF/ 6-31G level with those of Alonso et al. experimentally observed conformers of isolated α- and β-D-glucose, which used Fourier transform microwave spectroscopy coupled with laser ablation.64 In their experiments, four conformers were observed for α-D-glucose, which were labeled as G−g+/cc/t*, G+g−/cc/t*, Tg+/cc/t, and G−g+/cl/g−, and three conformers were observed for β-D-glucose, which were labeled as G−g+/cc/t*, G+g−/cc/t*, and Tg+/cc/t. These labels indicate plausible conformations of the hydroxymethyl groups (Figure 15): The

Figure 15. Conformer notations for hydroxylmethyl groups. (a) Newman projections of the plausible conformations of the hydroxylmethyl group around C5−C6. (b) Newman projections of the plausible conformations of the hydroxylmethyl group around C6− O6. (c) Orientation of the hydroxyl groups. (d) Newman projections of the plausible orientation of the anomer hydroxyl group hydrogen atom.

symbol G+, G−, or T designates the torsion angle O6−C6− C5−O5 of about 60°, −60°, or 180°, respectively (Figure 15a). The following symbol g+, g−, or t designates the torsion angle H6−O6−C6−C5 in the same way as shown in Figure 15b. The third symbol between slashes stands for the orientation of the hydroxyl groups: cc (counter clockwise) or cl (clockwise) (Figure 15c). The last symbol g+, g−, or t is regarding the torsion angle C2−C1−O1−H1 (Figure 15d). The most abundant conformer types observed in the experiments are marked with an asterisk. We examined if these experimentally obtained conformers were reproduced by the PES-based L


Article


Figure 16. Geometries of the five conformers of α-D-glucose in low energy giving less than 10 kJ/mol relative to 1amin obtained by the PES-based search at RHF/6-31G with LADD = 5. The intramolecular O···H distances shorter than 3.0 Å are shown in Å. In the parentheses, the energies are given relative to 1amin from the left side at the RHF/6-31G//RHF/6-31G, B3LYP/6-31G//RHF/6-31G, and MP2/6-311++(d,p)//RHF/6-31G levels (footnotes b, c, and e in Table 1). The notation below it designates the conformation for the hydroxylmethyl groups (Figure 15). The conformational notations with star mark(s) are ones that were also experimentally observed by Alonso et al. using Fourier transform microwave spectroscopy coupled with laser ablation of crystalline α-D-glucose;64 double stars are given to the most abundant conformers in the experiments.

Table 1. Relative Energies Corresponding to the Conformers in Low Potential Energy Obtained by the PES-based Search at RHF/6-31G with LADD = 5 relative energya [kJ/mol] entry no.

structure no.

HFb

B3LYP-singlec

B3LYP-optd

MP2-singlee

MP2-optf

conformation typeg

exp.h obsd: ×

1 2 3 4 5 6 7 8

1amin 1a2 1a3 1a4 1a5 2amin 2a2 2a3

0 2.0 3.1 4.0 6.0 0 0.1 2.1

0 −4.0 0.7 1.5 5.7 0 −3.4 1.8

0 −3.7 1.3 2.2 5.9 0 −3.7 1.9

0 −3.0 −4.4 −3.7 −1.9 0 2.1 0.9

0 −3.8 −5.6 −5.0 −2.2 0 2.3 0.7

Tt/cl/g− Tg+/cc/t G−g+/cc/t G+g-/cc/t G−g+/cl/g− G−g+/cc/t Tg+/cc/t G+g−/cc/t

× ×* ×* × ×* × ×*

The energies are given relative to entries 1 (for the α-anomer) or 6 (for the β-anomer). bRHF/6-31G//RHF/6-31G. cB3LYP/6-31G//RHF/631G. dB3LYP/6-31G//B3LYP/6-31G. eMP2/6-311++G(d,p)//RHF/6-31G. fMP2/6-311++G(d,p)// MP2/6-311++G(d,p). gConformational notation for hydroxylmethyl groups shown in Figure 15. hJ. L. Alonso et al. Chem. Sci. 2014, 5, 515−522. *: the most abundant conformation types in the experiments. a

Figure 17. Geometries of the three conformers of β-D-glucose in low energy giving less than 10 kJ/mol relative to 2amin obtained by the PES-based search at RHF/6-31G with LADD = 5. The intramolecular O···H distances shorter than 3.0 Å are shown. In the parentheses, the energies are given relative to 2amin from the left side at the RHF/6-31G//RHF/6-31G, B3LYP/6-31G//RHF/6-31G, and MP2/6-311++(d,p)//RHF/6-31G levels (footnotes b, c, and e in Table 1). The notation below it shows the conformation for the hydroxylmethyl groups (Figure 15). The conformational notations with star mark(s) are ones that were also experimentally observed by Alonso et al. using Fourier transform microwave spectroscopy coupled with laser ablation of crystalline β-D-glucose;64 double stars are given to the most abundant conformers in the experiments. M


Article

Journal of Chemical Theory and Computation same way. At the B3LYP/6-31G level, the experimentally least conformer (2a2) was estimated to be the minimum (2a2 → 2amin → 2a3), but excellent agreement was found at the MP2/ 6-311++G(d,p) level (2amin → 2a3 → 2a2). After the geometry optimization, only very slight changes were found in energies and geometries that hold the same conformation types in any level (Table 1 footnotes d and f, respectively). These results demonstrate that the PES-based conformational search at RHF/6-31G reproduced well the geometries observed in experiments. The RHF/6-31G geometries in low potential energy are good agreement with the B3LYP/6-31G and MP2/6-311++G(d,p) geometries. The RHF/6-31G energies are valid enough to use for qualitative discussions. For quantitative discussions, energy correction using a better QM method and a larger basis set after a search will be necessary. The QM method and the basis set used in the PESbased conformational search should be determined by considering the problem to solve, the size of the molecular system, and the computer resource available for the calculations. The verification results suggest that conducting a conformational search with a moderate QM level followed by an energy correction with a higher QM level can be a valid option for an efficient search.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS H.S. is very grateful to Professor Dennis M. Whitfield (National Research Council, Canada) for kindly providing us their program to calculate the CP coordinates. H.S. is also very thankful to Professor Jürg Hutter (University of Zurich, Switzerland) for valuable discussions. The project is partially supported by the Grant from the Data Centric Science Research Commons Project of the Research Organization of Information and Systems (ROIS), Japan. H.S. and T.U. were supported by a Grant-in-Aid for Challenging Exploratory Research (Grant No. 25540017) from the Japan Society for the Promotion of Science, and T.U. and K.N. are supported by CREST, Japan Science and Technology Agency (JST).

■

7. CONCLUSIONS This paper presented our new approach for a PES-based conformational analysis based on SHS-ADDF and the RMapViewer software we have developed. We demonstrated its applications to deducing conformational r-maps of α- and βD-glucoses. The PES-based conformational search has been carried out with large ADDF while restraining bond lengths and structures of high free energy in the exploration with SHSADDF. It automatically deduced conformational transition networks, called conformational r-maps, which are global conformational reaction route maps composed of EQ conformers connected by reaction routes via TS conformers. Postcalculation data analysis of the conformational r-maps with RMapViewer revealed MEPs between the global minima of 4C1 and 1C4 conformations. The conformers found in low potential energy by the PES-based search are in good agreement with the experiments by Alonso et al.64 Although this paper focused our discussions on MEPs in the conformational r-maps, our analysis using RMapViewer has also revealed that the obtained r-maps contain a larger number of reaction routes with closer TSmax values to those of MEPs. If those reaction routes are bundled, a conformational reaction flow map becomes derivable, which would look like a landscape of reaction streams. Furthermore, an r-map may contain new reaction routes, which have not been discussed yet but would likely to have significant impact in chemical or biochemical reactions. The bond-breaking channels in low energy barriers obtained by SHS-ADDF with the restrictions would include useful transformations in the synthetic design of glucopyranosides or in the analysis of bioprocess of glucose. R-maps, together with RMapViewer, will provide a vast amount of information that would be helpful in analyzing and in designing of molecules and chemical reaction routes.

■

All of the input files to GRRM and conformer types plotted against energy obtained with LADD = 10 (PDF)

REFERENCES

(1) Leach, A. R. Conformational Analysis. Molecular Modelling. Principles and Applications, 2nd ed.; Pearson Education Limited: Edinburgh Gate Harlow Essex CM20 2JE, England, 2001; pp 457− 508. (2) Jensen, F. Conformational Sampling and the Global Minimum Problem. Introduction to Computational Chemistry, 2nd ed.; John Wiley & Sons Ltd.: The Atrium, Southern Gate, Chichester, West Sussex PO19 8SQ, England, 2007; pp 409−415. (3) For example, if the systematic rotation is done even only for 10 rotatable bonds in steps of 30°, 1210 (= 61,917,364,224) conformers are explored. (4) Chang, G.; Guida, W. C.; Still, W. C. An Internal Coordinate Monte Carlo Method for Searching Conformational Space. J. Am. Chem. Soc. 1989, 111, 4379−4386. (5) Chandrasekhar, J.; Saunders, M.; Jorgensen, W. Efficient Exploration of Conformational Space Using the Stochastic Search Method: Application to β-Peptide Oligomers. J. Comput. Chem. 2001, 22, 1646−1654. (6) Huber, T.; Torda, A. E.; van Gunsteren, W. F. Local elevation: A method for improving the searching properties of molecular dynamics simulation. J. Comput.-Aided Mol. Des. 1994, 8, 695−708. (7) Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562−12566. (8) Henkelman, G.; Jónsson, H. Improved Tangent Estimate in The Nudged Elastic Band Method for Finding Minimum Energy Paths and Saddle Points. J. Chem. Phys. 2000, 113, 9978−9985. (9) Henkelman, G.; Uberuaga, B. P.; Jónsson, H. A Climbing Image Nudged Elastic Band Method for Finding Saddle Points and Minimum Energy Paths. J. Chem. Phys. 2000, 113, 9901−9904. (10) Zarkevich, N. A.; Johnson, D. D. Nudged-elastic Band Method with Two Climbing Images: Finding Transition States in Complex Energy Landscapes. J. Chem. Phys. 2015, 142, 024106. (11) Krivov, S. V.; Karplus, M. Hidden Complexity of Free Energy Surfaces for Peptide (Protein) Folding. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 14766−14770. (12) Noé, F.; Fischer, S. Transition Networks for Modeling The Kinetics of Conformational Change in Macromolecules. Curr. Opin. Struct. Biol. 2008, 18, 154−162. (13) Trygubenko, S. A.; Wales, D. J. A Doubly Nudged Elastic Band Method for Finding Transition States. J. Chem. Phys. 2004, 120, 2082. (14) Carr, J. M.; Trygubenko, S. A.; Wales, D. J. Finding Pathways Between Distant Local Minima. J. Chem. Phys. 2005, 122, 234903.

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b00439. N


Article

Journal of Chemical Theory and Computation (15) Wales, D. J. Energy Landscapes: Some New Horizons. Curr. Opin. Struct. Biol. 2010, 20, 3−10. (16) Grouleff, J.; Jensen, F. Searching Peptide Conformational Space. J. Chem. Theory Comput. 2011, 7, 1783−1790. (17) Halgren, T. A.; Lipscomb, W. N. The Synchronous-transit Method for Determining Reaction Pathways and Locating Molecular Transition States. Chem. Phys. Lett. 1977, 49, 225−232. (18) Bell, S.; Crighton, J. S. Locating Transition States. J. Chem. Phys. 1984, 80, 2464−2475. (19) Ohno, K.; Maeda, S. A Scaled Hypersphere Search Method for the Topography of Reaction Pathways on the Potential Energy Surface. Chem. Phys. Lett. 2004, 384, 277−282. (20) Maeda, S.; Ohno, K. Ab initio Studies on Synthetic Routes of Glycine from Simple Molecules via Ammonolysis of Acetolactone: Applications of the Scaled Hypersphere Search Method. Chem. Lett. 2004, 33, 1372−1373. (21) Maeda, S.; Ohno, K. No Activation Barrier Synthetic Route of Glycine from Simple Molecules (NH 3 , CH 2 , and CO 2 via Carboxylation of Ammonium Ylide: A Theoretical Study by the Scaled Hypersphere Search Method. Chem. Phys. Lett. 2004, 398, 240−244. (22) Maeda, S.; Ohno, K. A New Approach for Finding a Transition State Connecting a Reactant and a Product without Initial Guess: Applications of the Scaled Hypersphere Search Method to Isomerization Reactions of HCN, (H2O)2, and Alanine Dipeptide. Chem. Phys. Lett. 2005, 404, 95−99. (23) Kolossváry, I.; Guida, W. C. Low Mode Search. Low Mode Search. An Efficient, Automatic Computational Method for Conformational Analysis: Application to Cyclic and Acyclic alkanes and Cyclic Peptides. J. Am. Chem. Soc. 1996, 118, 5011−5019. (24) Kishimoto, N.; Harayama, M.; Ohno, K. An Automated Efficient Conformation Search of L-serine by the Scaled Hypersphere Search Method. Chem. Phys. Lett. 2016, 652, 209−215. (25) Maeda, S.; Ohno, K. Global Mapping of Equilibrium and Transition Structures on Potential Energy Surfaces by the Scaled Hypersphere Search Method: Application to ab initio Surfaces of Formaldehyde and Propyne Molecules. J. Phys. Chem. A 2005, 109, 5742−5753. (26) Yang, X.; Maeda, S.; Ohno, K. Global Investigation on Potential Energy Surface of CH3CN: Application of the Scaled Hypersphere Search Method. J. Phys. Chem. A 2005, 109, 7319−7328. (27) Maeda, S.; Watanabe, Y.; Ohno, K. A Scaled Hypersphere Interpolation Technique for Efficient Construction of Multidimensional Potential Energy Surfaces. Chem. Phys. Lett. 2005, 414, 265− 270. (28) Maeda, S.; Ohno, K.; Morokuma, K. Systematic Exploration of the Mechanism of Chemical Reactions: Global Reaction Route Mapping (GRRM) Strategy by the ADDF and AFIR Methods. Phys. Chem. Chem. Phys. 2013, 15, 3683−3701. (29) Ohno, K. Study of Potential Energy Surface towards Global Reaction Route Mapping. Chem. Rec., published online on 5 April 2016, DOI: 10.1002/tcr.201500284. ” is Japanese, which is a prefix for (30) The term “Maizo buried precious items. This term is generally used for buried precious items, such as cultural properties, fossils, treasure troves, and deposits of precious metals, oil, or coal. For example, “Maizo”-gold is a usage of this word. We used this word as part of the name of our project since it clearly expresses the concept towards finding out new chemistry buried in the QM-based big data space. (31) Satoh, H.; Oda, T.; Nakakoji, K.; Uno, T.; Iwata, S.; Ohno, K. Maizo”-chemistry Project: Towards Molecular- and Reaction Discovery from Quantum Mechanical Global Reaction Route Mappings. J. Comput. Chem., Jpn. 2015, 14, 77−79. (32) Ohno, K.; Satoh, H.; Iwamoto, T. A Prism Carbon Molecule C20. Chem. Lett. 2015, 44, 712−714. (33) Ohno, K.; Satoh, H.; Iwamoto, T. Prism-C2n Carbon Dimer, Trimer, and Nano-Sheets: A Quantum Chemical Study. Chem. Phys. Lett. 2015, 633, 120−125.

(34) Ohno, K.; Satoh, H.; Iwamoto, T.; Tokoyama, H.; Yamakado, H. Wavy Carbon: A New Series of Carbon Structures Explored by Quantum Chemical Calculations. Chem. Phys. Lett. 2015, 639, 178− 182. (35) Satoh, H.; Oda, T. RMapViewer, Version 4.0. RMapViewer is free to download from http://sourceforge.net/projects/rmapviewer/ (accessed August 19, 2016). The software is built on Pharo, an open source implementation of the programming language and environment Smalltalk and runs on Linux, Windows, and Mac OS. We have been preparing to release RMapViewer also as an open source code. (36) For example, an exploration of C6H6 isomers in the SCC-DFTB level with GRRM gave 7,000 isomers and 26,229 transition structures. See: Tokoyama, H.; Yamakado, H.; Maeda, S.; Ohno, K. Exploration of Isomers of Benzene by using the GRRM/SCC-DFTB Method. Chem. Lett. 2014, 43, 702−704. (37) Angyal, S. J. Conformational Analysis in Carbohydrate Chemistry I. Conformational Free Energies. The Conformations and α:β Ratios of Aldopyranoses in Aqueous Solution. Aust. J. Chem. 1968, 21, 2737−46. (38) Miljković, M. Conformational Analysis of Monosaccharides. Carbohydrate. Synthesis, Mechanisms, and Stereoelectronic Effects, 1st ed.; Springer Science+Business Media, LLC: 233 Spring Street, New York, NY 10013, USA, 2009; pp 27−56. (39) Satoh, H.; Manabe, S. Design of Chemical Glycosyl Donors: does changing ring conformation influence selectivity/reactivity? Chem. Soc. Rev. 2013, 42, 4297−4309. (40) Satoh, H.; Nukada, T. Computational Chemistry on Chemical Glycosylations. Trends Glycosci. Glycotechnol. 2014, 26, 11−27. (41) Bohé, L.; Crich, D. A propos of glycosyl cations and the mechanism of chemical glycosylation; the current state of the art. Carbohydr. Res. 2015, 403, 48−59. (42) Whitfield, D. M. In a Glycosylation Reaction How Does a Hydroxylic Nucleophile Find The Activated Anomeric Carbon? Carbohydr. Res. 2015, 403, 69−89. (43) Kirby, A. J. Effects on conformation. Stereoelectronic Effects; Oxford University Press: Great Clarendon Street, Oxford ox2 6DP, UK, 1996; pp 13−23. (44) Miljković, M. Anomeric Effect. Carbohydrate. Synthesis, Mechanisms, and Stereoelectronic Effects, 1st ed.; Springer Science +Business Media, LLC: 233 Spring Street, New York, NY 10013, USA, 2009; pp 57−93. (45) Melberg, S.; Rasmussen, K.; Scordamaglia, R.; Tosi, C. The Non-bonded Interactions in D-Glucose and β-Maltose: An ab initio Study of Conformations Produced by Empirical Force-Field Calculations. Carbohydr. Res. 1979, 76, 23−37. (46) Cramer, C. J.; Truhlar, D. G. Quantum Chemical Conformational Analysis of Glucose in Aqueous Solution. J. Am. Chem. Soc. 1993, 115, 5745−5753. (47) Cramer, C. J.; Truhlar, D. G.; French, A. D. Exo-anomeric Effects on Energies and Geometries of Different Conformations of Glucose and Related Systems in The Gas Phase and Aqueous Solution. Carbohydr. Res. 1997, 298, 1−14. (48) Ionescu, A. R.; Bérces, A.; Zgierski, M. Z.; Whitfield, D. M.; Nukada, T. Conformational Pathways of Saturated Six-Membered Rings. A Static and Dynamical Density Functional Study. J. Phys. Chem. A 2005, 109, 8096−8105. (49) Biarnés, X.; Ardèvol, A.; Planas, A.; Rovira, C.; Laio, A.; Parrinello, M. The Conformational Free Energy Landscape of β-DGlucopyranose. Implications for Substrate Preactivation in β-Glucoside Hydrolases. J. Am. Chem. Soc. 2007, 129, 10686−10693. (50) Schnupf, U.; Willett, J. L.; Momany, F. DFTMD Studies of Glucose and Epimers: Anomeric Ratios, Rotamer Populations, and Hydration Energies. Carbohydr. Res. 2010, 345, 503−511. (51) Barnett, C. B.; Wilkinson, K. A.; Naidoo, K. J. Pyranose Ring Transition State Is Derived from Cellobiohydrolase I Induced Conformational Stability and Glycosidic Bond Polarization. J. Am. Chem. Soc. 2010, 132, 12800−12803. O


Article

Journal of Chemical Theory and Computation (52) Andrade, R. R.; da Silva, C. O. Specific Rotation as a Property to Validate Monosaccharide Conformations. Carbohydr. Res. 2012, 350, 62−67. (53) Biarnés, X.; Ardèvol, A.; Iglesias-Fernández, J.; Planas, A.; Rovira, C. Catalytic Itinerary in 1,2−1,4−β−Glucanase Unraveled by QM/MM Metadynamics. Charge Is Not Yet Fully Developed at the Oxacarbenium Ion-like Transition State. J. Am. Chem. Soc. 2011, 133, 20301−20309 See also ref 9 therein.. (54) Davies, G. J.; Planas, A.; Rovira, C. Conformational Analyses of the Reaction Coordinate of Glycosidases. Acc. Chem. Res. 2012, 45, 308−316 See also refs 22, 23, and 51 therein.. (55) Mayes, H. B.; Broadbelt, L. J.; Beckham, G. T. How Sugars Pucker: Electronic Structure Calculations Map the Kinetic Landscape of Five Biologically Paramount Monosaccharides and Their Implications for Enzymatic Catalysis. J. Am. Chem. Soc. 2014, 136, 1008−1022. (56) Maeda, S.; Osada, Y.; Morokuma, K.; Ohno, K. GRRM, GRRM11. To use the GRRM program, directly contact K.O., the representative of the developers. See the detailed information about the distribution of the GRRM program package at http://grrm.chem. tohoku.ac.jp/GRRM/index_e.html (accessed August 19, 2016). Information on the newest version will be appearing. (57) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, Revision D.01; Gaussian, Inc.: Wallingford, CT, 2009. (58) Maeda, S.; Ohno, K. Structures of Water Octamers (H2O)8: Exploration on Ab Initio Potential Energy Surfaces by the Scaled Hypersphere Search Method. J. Phys. Chem. A 2007, 111, 4527−4534. (59) Satoh, H.; Koshino, H.; Funatsu, K.; Nakata, T. A Novel Canonical Coding Method for Representation of Three-dimensional Structures. J. Chem. Inf. Comput. Sci. 2000, 40, 622−630. (60) Satoh, H.; Koshino, H.; Funatsu, K.; Nakata, T. Representation of Configurations by CAST Coding Method. J. Chem. Inf. Comput. Sci. 2001, 41, 1106−1112. (61) IUPAC−IUB Joint Commission on Biochemical Nomenclature (JCBN). Conformational Nomenclature for Five and Six-membered Ring Forms of Monosaccharides and Their Derivatives. Eur. J. Biochem. 1980, 111, 295−298; Arch. Biochem. Biophys. 1981, 207, 469−472; Pure Appl. Chem. 1981, 53, 1901−1905. (62) Cremer, D.; Pople, J. A. A General Definition of Ring Puckering Coordinates. J. Am. Chem. Soc. 1975, 97, 1354−1358. (63) Characterization of Six-membered Ring Conformations: http:// 6ring.bio.nrc.ca (accessed August 19, 2016). (64) Alonso, J. L.; Lozoya, M. A.; Peña, I.; López, J. C.; Cabezas, C.; Mata, S.; Blanco, S. The Conformational Behaviour of free D-glucose − at last. Chem. Sci. 2014, 5, 515−522.

P


Potential Energy Surface-Based Automatic Deduction of

Recommend Documents