Low-Barrier Hydrogen Bond in Fujikurin A-D: A Computational Study.

The compounds Fujikurin A, B, and D, recently isolated from Fusarium fujikuroi, possess intramolecular low-barrier hydrogen bonds (LBHBs), which are hydrogen bonds with a very low-energy barrier for proton transfer. The isolated compounds have a hydrogen-bonded proton that appears to rapidly switch between two equilibrium states via a transition state (TS). To understand the characteristics of these intramolecular LBHBs in detail, we performed path integral molecular dynamics (PIMD) simulations, which can consider nuclear quantum effects (NQEs) under a finite temperature. The PIMD simulations predicted that the NQE completely washed out the energy barrier for the proton transfer reaction. Consequently, a single-well shape emerged in the results, along with the effective free-energy potential surface for the hydrogen-bonded proton distribution. Thus, we conclude that the hydrogen-bonded proton in Fujikurin does not in fact transfer between two equilibrium structures but widely delocalizes around the global minimum structure involving the TS region.

Introduction

Fungi of the genus Fusarium, including Fusarium fujikuroi, produce harmful substances that include mycotoxins. They also induce crop damage with gibberellins that cause the rice disease bakanae. Fusarium fujikuroi is known to produce various secondary metabolites in addition to gibberellins.[1,2] These metabolites have been studied extensively over time.[1−8] Three new compounds, Fujikurin A, B, and D (Figure ), were recently identified by Bargen et al. These compounds are a secondary metabolite produced by polyketide synthases (PKSs).[6]

Figure 1

Chemical structures of Fujikurin A–D.

Bargen et al. focused on the main product, Fujikurin D, analyzing its structure experimentally (mass spectrometry and NMR) and theoretically (density functional theory, DFT), and determined that it had an intramolecular hydrogen bond. Furthermore, they proposed that the hydrogen-bonded proton may rapidly shift back and forth between two equilibrium structures, EQ(1) and EQ(2) (Figure ). Moreover, Fujikurin A and B were also found to have intramolecular hydrogen bonds that have the same characteristics as those of Fujikurin D.

Figure 2

Equilibrium and TS structures of Fujikurin D.

However, most of the widely used theoretical calculation methods based on the Born–Oppenheimer approximation do not consider nuclear motion (nuclear quantum effects (NQEs)). The fixed nuclear approach based on the Born–Oppenheimer approximation works well in most situations due to the fact that nuclei are considerably heavier than electrons. However, to analyze systems containing strong hydrogen bonds, including low-barrier hydrogen bonds (LBHBs), where the energy barrier for proton transfer is very low, NQEs are highly useful, enabling the qualitative reproduction of structures.[9−12] Strongly hydrogen-bonded systems have been studied with many helpful methods, including those that involve NQEs. For example, Litman et al. studied porphycene, which has strong intramolecular hydrogen bonds, and reported the temperature dependence of the double hydrogen transfer dynamics, which can be well-understood only if NQEs are taken into account.[13] Ogata et al. investigated asymmetric LBHBs in protonated lysine (LysH+), analyzing the NQEs using the path integral molecular dynamics (PIMD) method.[11] They determined that the hydrogen-bonded proton is not localized on the proton donor or acceptor heavy atoms but is instead delocalized between the heavy atoms from the PIMD simulation.

Thus, the consideration of NQEs is important to correctly understand LBHBs. However, Fujikurin D was analyzed via conventional DFT methods in Bargen’s work, which means that NQEs were not included. Many theoretical methods have been proposed to include NQEs. For example, we have developed multicomponent molecular orbital (MC_MO) methods[9,10,14−16] that go beyond the Born–Oppenheimer approximation to extend the concept of molecular orbitals for electrons to light nuclei. We used MC_MO methods to successfully evaluate the effects of NQEs and H/D isotopes in different hydrogen-bonded systems and proton transfer reactions. PIMD can also consider NQEs statistically within the Born–Oppenheimer approximation framework.[11−13,17−24] The PIMD approach considers NQEs as well as thermal effects. Many PIMD studies have been used to understand NQEs in the context of hydrogen-bonded systems, including LBHBs. In this study, we concentrated on the distribution of hydrogen-bonded protons within Fujikurin A–D. Hence, we carried out PIMD simulations on Fujikurin A–D. We also incorporated typical DFT calculations and conventional MD simulations to elucidate the NQEs and thermal effects in Fujikurin A–D.

Computational Methods

Nuclei are treated as quantum particles in PIMD simulations. In this study, each nucleus was expanded using 16 beads. PIMD simulations were then performed for 90 000 steps, after a thermal equilibration of 10 000 steps using a time step size of 0.1 fs. Furthermore, classical MD (CLMD) simulations, in which each nucleus was represented by one bead, were carried out for 900 000 steps following a thermal equilibration of 100 000 steps using the same time step size for comparison. We performed both CLMD and PIMD simulations in the canonical ensemble (NVT) using a Nosé–Hoover chain[31] at 300 K.

The position of the hydrogen-bonded proton was clarified by defining the proton transfer coordinate (δOH) aswhere RO1H and RHO2 are the O1–H covalent bond length and H···O2 hydrogen bond distance, respectively. The notation is illustrated in Figure . The value of δOH = 0 implies that the proton is at the center of the oxygen atoms, while a positive or negative δOH value means that the proton is located at a position closer to one of the oxygen atoms.

Figure 3

Hydrogen bond structure of Fujikurin D and definition of atomic labels and coordinates.

The relative effective free energy (ΔF) was calculated from the obtained distributions using the following equation:where kB, T, and P(δOH) are the Boltzmann constant, temperature, and probability distribution, respectively, as a function of δOH. P(δOHmax)) is the largest distribution of P(δOH) in the histogram, and the free energy in this region is set to 0.0 kcal/mol.

Results and Discussion

Level of Electronic Structure Calculation

The level of electronic structure calculations for the PIMD simulations was based on the relative energy of each stationary point structure of the main product, Fujikurin D. Table shows the relative energies obtained by DFT (CAM-B3LYP),[25] MP2, and semiempirical PM series methods.[26−30]

Table 1

Relative Energies (kcal/mol) of Stationary Point Structures of Intramolecular Proton Transfer Reaction in Fujikurin D Obtained by Several Methods

method	EQ(1)	TS	EQ(2)
CAM-B3LYP/6-311+G**	0.0	1.8	1.1
MP2/aug-cc-pVDZ//CAM-B3LYP/6-311+G**	0.0	1.9	0.9
PM6	0.0
PM6-DH2	0.0
PM6-DH+	0.0	3.2	2.0
PM6-D3H4	0.0	3.0	1.9
PM7	0.0	2.7	2.5

The MP2/aug-cc-pVDZ//CAM-B3LYP/6-311+G** energy was considered to be a highly accurate reference. The CAM-B3LYP/6-311+G** energies were similar to the MP2/6-311+G** energies. Hence, the CAM-B3LYP/6-311+G** method describes the potential energy surface (PES) for the proton-transfer reaction in Fujikurin D accurately. It was also determined that EQ(1) has a lower energy than that of EQ(2). Therefore, EQ(1) is the global minimum structure, and EQ(2) is the metastable equilibrium structure. PIMD simulations require very large numbers of atomic force calculations; therefore, selecting a method that can reproduce highly accurate results with a low computational cost is very important. We focused on a semiempirical PM series method. In Table , only the equilibrium state EQ(1) exists on the PM6 and PM6-DH2 PESs, while the proton-transfer equilibrium state EQ(2) and the TS exist on the PM6-DH+, PM6-D3H4, and PM7 PESs. In particular, the PM6-D3H4 method reproduced the relative energies obtained by the MP2/aug-cc-pVDZ//CAM-B3LYP/6-311+G** method.

Next, we confirmed the validity of the PM6-D3H4 method for Fujikurins A and B. The relative energies of Fujikurin A and B are listed in Table .

Table 2

Relative Energies (kcal/mol) of Stationary Point Structures of Intramolecular Proton-Transfer Reaction in Fujikurin A and B Obtained by Several Methods

	Fujikurin A			Fujikurin B
method	EQ(1)	TS	EQ(2)	EQ(1)	TS	EQ(2)
CAM-B3LYP/6-311+G**	0.0	1.4	0.2	0.0	2.3	1.7
MP2/aug-cc-pVDZ//CAM-B3LYP/6-311+G**	0.0	1.4	–0.2	0.0	2.2	1.4
PM6-D3H4	0.0	2.8	1.7	0.0	3.7	2.9

Similar to the case of Fujikurin D, the CAM-B3LYP/6-311+G** relative energies are like those calculated by MP2/aug-cc-pVDZ//CAM-B3LYP/6-311+G**, indicating that the CAM-B3LYP-optimized geometries were adequate for Fujikurin A and B. The MP2/aug-cc-pVDZ//CAM-B3LYP/6-311+G** activation barriers of Fujikurin A and Fujikurin B were 0.5 kcal/mol lower and 0.3 kcal/mol higher than that of Fujikurin D, respectively. Moreover, the PM6-D3H4 method estimated the activation barriers of Fujikurin A and Fujikurin B as 0.2 kcal/mol lower and 0.7 kcal/mol higher than those of Fujikurin D, respectively. This implies that the PM6-D3H4 method is also adequate for analyzing the intramolecular proton-transfer reaction in Fujikurin A and B. Finally, we used the PM6-D3H4 method to calculate atomic forces.

Structural Features of Fujikurin D

For Fujikurin D, Figure shows the one-dimensional distributions of RO1H, RHO2, and RO1O2, while Table has the average values of RO1H, RHO2, and RO1O2 obtained by CLMD and PIMD simulations. The values shown in parentheses are the statistical errors of the average bond lengths calculated using the block average method.[32]Table lists the optimized RO1H, RHO2, and RO1O2 values in the equilibrium structure.

Figure 4

One-dimensional distributions of (a) RO1H, (b) RHO2, and (c) RO1O2 obtained by CLMD and PIMD simulations. The dashed lines represent the values of two equilibrium structures.

Table 3

Average Values of RO1H (Å), RHO2 (Å), and RO1O2 (Å) Obtained by CLMD and PIMD Simulations and the Interatomic Distances in Equilibrium Structures EQ(1) and EQ(2)

	RO1H	RHO2	RO1O2
CLMD	1.103(7)	1.607(8)	2.525(2)
PIMD	1.112(4)	1.601(10)	2.512(6)
EQ(1)	1.091	1.582	2.504
EQ(2)	1.420	1.174	2.440

In a comparison of the CLMD and PIMD results for the one-dimensional distributions of RO1H and RHO2 (Figure a,b), the distributions generated by PIMD simulations are much wider than those produced by CLMD simulations owing to quantum fluctuations. The results of the CLMD simulations showed two peaks. Because EQ(1) is more stable than EQ(2), the peak around EQ(1) is much larger. However, the distributions for the CLMD and PIMD simulations for RO1O2 are similar. Thus, the impact of the NQEs on RO1O2 was smaller than that on RO1H.

We focused on the covalent RO1H average bond lengths in Table and found that the relationship of the RO1H values is as follows:

The average value obtained in the CLMD simulation was greater than the ROH in EQ(1). The CLMD simulation includes thermal effects only; therefore, the elongation of the covalent RO1H bond length is due to thermal motion. The PIMD simulation, which includes both thermal effects and NQEs, provided a greater RO1H average bond length than that of the CLMD simulation. Thus, the NQEs stretch the covalent bond length. These results are often observed in other hydrogen-bonded systems.[11,12]

However, the average values of RHO2 and RO1O2 in Table show the following relationship, which is different from that for the covalent ROH lengths:

The bond lengths obtained from the CLMD and PIMD simulations were greater than those in EQ(1) because of the thermal effect, as in the case of RO1H. To explain why the RHO2 and RO1O2 values obtained in the PIMD simulations were smaller than those obtained in the CLMD simulation, we focused on the Mulliken charges of oxygen and hydrogen atoms (Figure ).

Figure 5

Mulliken charges of the oxygen and the hydrogen atoms obtained in (a) CLMD and (b) PIMD simulations.

The Mulliken charge of H was found to be +0.426 in the CLMD result, whereas it was +0.430 from PIMD, which means that the hydrogen-bonded proton in the PIMD simulation was more positively charged than that in the CLMD simulation. The finite spread of the nuclei in the PIMD simulation induced the NQEs, which affected the surrounding electron distributions. We can find that the proton acceptor, O2, is more negative than O1. In addition, O2 in the PIMD simulation is less negative than that in the CLMD one, unlike O1. Since the proton donor and acceptor atoms swap their roles by the proton transfer, O1 and O2 become more and less negative around the TS region. The hydrogen-bonded proton can access the TS region more easily in the PIMD simulation rather than in the CLMD one, as described below. Thus, O2 seems to be less negative in the PIMD simulation than in the CLMD one. We should focus on the sum of the negative charges of O1 and O2 for the strength of the interaction. The greater negative charge is found in the PIMD simulation than in the CLMD. The more positively charged hydrogen-bonded proton interacts more strongly with the negatively charged oxygen atoms in the hydrogen bond moiety. Thus, the distance between the oxygen atoms in the PIMD result was smaller than that in the CLMD result owing to the strong hydrogen bond interaction.

Relative Free-Energy Potential Surface for Proton-Transfer Coordinate

To understand the hydrogen-bonded proton distribution, we focused on the one-dimensional relative free-energy potential surface for the proton-transfer coordinate δOH. Figure shows the one-dimensional relative effective free-energy potential surfaces from the CLMD and PIMD simulations.

Figure 6

One-dimensional relative effective free-energy landscapes of δOH obtained by CLMD and PIMD simulations. The dashed lines represent two equilibrium structures.

First, we focused on the CLMD results. The relative effective free-energy potential obtained by the CLMD simulation is a double-well potential with two minima. The positions of the minima are in good agreement with the δOH values of the equilibrium structures shown by the black dashed lines in Figure . Thus, it is difficult to overcome the activation energy barrier using thermal effects alone. The barrier height of the free-energy potential in the CLMD result was estimated as 3.6 kcal/mol, which is slightly larger than that from the static electronic structure calculation, 3.0 kcal/mol (Table ). This trend has also been reported for other LBHB systems.[11,17] However, the energy barrier completely disappeared in the PIMD results, and the shape of the free-energy potential was a single well. The PIMD simulation took the thermal effect and NQEs into account; thus, the energy barrier of the proton transfer was washed out by the NQEs. Similar NQE influences on the proton distribution around the TS structure were reported for LBHB in LysH+ by Ogata et al.[11] In particular, the NQEs on proton transfer might be the primary factor because the NQEs of light nuclei were more prominent than those of heavier nuclei. It should be noted here that, in a symmetrical hydrogen-bonded system such as a maleate anion,[12] the potential minimum should be at δOH = 0 Å. However, in Fujikurin D, the minimum potential is located at approximately δOH = −0.5 Å, which corresponds to the global minimum structure EQ(1), because Fujikurin D has an asymmetric PES in the proton-transfer coordinate.

The two-dimensional free-energy landscapes of δOH and ROO are shown in Figure to aid a better understanding of the proton distribution. There are two minima located around the equilibrium structures in the free-energy landscape simulated using the CLMD method (Figure a). Thus, the equilibrium structure analogues, in the which transferring proton is localized to either oxygen atom, are frequently observed as expected from the one-dimensional free-energy landscape (Figure ). In addition, the energy maximum exists near the TS structure. Therefore, the distribution near the TS region is sparse compared to those around the equilibrium structures. Moreover, the low-energy regions are distributed along the IRC pathways. Thus, the hydrogen-bonded proton is mainly distributed around the equilibrium structures, and proton transfer may occur between the two equilibrium structures. Thus, the predictions of the CLMD are consistent with those by Bargen.

Figure 7

Two-dimensional free-energy landscapes of δOH and ROO obtained with (a) CLMD and (b) PIMD. Triangles and circles represent the equilibrium and the TS structures, respectively, and dashed lines are IRC pathways. Free energy is in units of kcal/mol.

In contrast, the PIMD simulation results are qualitatively different from the CLMD results. As seen in Figure b, the energy barrier for the proton transfer was completely overcome by the NQEs. The PIMD results show a broad single-well potential centered at the global minimum (δOH = −0.5 Å) and the corresponding one-dimensional free-energy potential (Figure ). The regions around the TS and metastable structures are also included in this broad single-well potential. Unlike in the CLMD predictions, the proton distribution extends to regions far from the IRC. This corner-cutting phenomenon is induced by the NQEs and is seen in several systems, such as N2H7+.[18] Our PIMD simulation shows a broad single-well potential centered on the global minimum structure for Fujikurin D. In other words, in the PIMD predictions, the proton overcomes the proton-transfer barrier freely and is not localized around either oxygen atoms, in contrast to Bargen’s prediction of fast transfer between the two equilibrium structures.

Proton Coordinate

We defined a coordinate system to analyze how the hydrogen-bonded protons are distributed, as shown in Figure . The origin is the center of mass of two oxygen atoms, and the xy-plane is defined to be in the O1–O2–C plane. Proton-distribution landscapes were obtained in three dimensions from both CLMD (Figure a) and PIMD (Figure b). Simulations in this coordinate system are shown in Figure .

Figure 8

Three-dimensional distribution of hydrogen-bonded proton in Fujikurin D. The origin is defined as the center of mass of the O atoms and the xy-plane as the O1–O2–C plane.

The CLMD results are shown in Figure a. Many structures are located around the two equilibrium structures, which were indicated by two peaks exhibited by the proton distribution in the xy- and xz-planes. In contrast, the PIMD result (Figure b) shows single broad distributions in every direction, depicting the activation barrier wash out caused by the NQEs, as mentioned above. While there are pronounced differences in the x- and y-directions when PIMD is compared to CLMD, the changes in the z-axis direction are not significant. To illustrate the changes more clearly, we show the difference distribution landscapes, which are the differences between the proton distributions obtained in the CLMD and PIMD simulations, in Figure .

Figure 9

Difference distribution landscapes of proton distributions obtained from CLMD and PIMD simulations of Fujikurin D. Red (blue) indicates regions where the distribution in the PIMD (CLMD) simulated results is denser than that in the CLMD (PIMD) simulated results.

In Figure , red indicates regions where the distributions in the PIMD results are denser than those in the CLMD results. This red area is spread in the xy-plane. Thus, distributions in the x- and y-directions increase because of the NQEs. However, in the yz- and zx-planes, the distribution in the z-direction did not increase much compared to the x- and y-directions. This means the NQEs did not considerably affect the protonic motion along the z-axis direction. The presence of the intramolecular hydrogen bond is considered to be a reason for the directional behavior of the NQE. The NQEs influence the electronic distribution, and as a result, the hydrogen bond interaction was enhanced in the PIMD simulations. Proton donor and acceptor atoms are always in the xy-plane in this environment. Hence, the hydrogen-bonded proton was attracted more strongly from the x- and y-directions than from the z-direction. This explains why the NQEs are pronounced in only the x- and y-directions.

Proton Distribution in Fujikurin A and B

Figure a,b shows the two-dimensional free-energy values, δOH and ROO, for Fujikurin A and B, respectively, to demonstrate how the distribution of the transferring proton is affected by the height of the activation barrier.

Figure 10

Two-dimensional free-energy landscapes δOH and ROO of Fujikurin A and B for (a) CLMD and (b) PIMD. Triangles and circles represent equilibrium and the TS structures, respectively, and dashed lines represent IRC pathways. Free energy is expressed in units of kcal/mol.

The CLMD results exhibit two minima, while the PIMD results show one minimum. This means the NQEs washed out the activation barrier in all Fujikurin compounds. Among the Fujikurin A–D compounds, Fujikurin A has the lowest activation barrier and largest low-energy region (represented by orange and yellow) in the two-dimensional free-energy results obtained by PIMD. NQEs extend the low-energy region around the TS, as well as near the metastable equilibrium structure.

It is worth noting that the PM6-D3H4 calculations overestimate the activation barrier slightly, when compared to those of the MP2/aug-cc-pVDZ//CAM-B3LYP/6-311+G** method. Hence, the distributions near the TS and metastable equilibrium structure may be denser on the more precise potential energy hypersurface.

Conclusion

In this study, we investigated the NQE and thermal effects on Fujikurin A–D intramolecular hydrogen bonds via PIMD simulation, which takes the thermal effect and the NQE into account.

For Fujikurin D, the CLMD simulations demonstrate that the hydrogen-bonded proton is rapidly switching between two equilibrium states, as predicted in a previous study by Bargen et al. In contrast, the PIMD simulation results are qualitatively different because NQEs are taken into consideration. The results imply a delocalized distribution on the single-well potential as a result of the wash out of the proton-transfer energy barrier. Furthermore, the NQE in Fujikurin D was found to be direction-dependent. The PIMD simulations demonstrated that the change in the PES of the proton-transfer reaction in Fujikurin D after considering NQEs could be observed even qualitatively.

Among Fujikurin A, B, and D, the widest proton distribution was found in Fujikurin A, which has the lowest activation barrier for proton transfer. In Fujikurin A, protons are distributed around the TS, as well as in the metastable equilibrium structure. Because the activation barriers were overestimated slightly by the PM6-D3H4 level of calculation, the proton distribution near the TS should actually be denser than calculated in this study. Our results can help elucidate the structures and properties of new secondary metabolites of Fujikurins.