2014 Volume 62 Issue 3 Pages 229-237
Compounds with a medium-sized flexible ring often show atropisomerism that is caused by the high-energy barriers between long-lived conformers that can be isolated and often have different biological properties to each other. In this study, the frequency of the transition between the two stable conformers, aS and aR, of thienotriazolodiazepine compounds with flexible 7-membered rings was estimated computationally by Monte Carlo (MC) simulations and validated experimentally by NMR experiments. To estimate the energy barriers for transitions as precisely as possible, the potential energy (PE) surfaces used in the MC simulations were calculated by molecular orbital (MO) methods. To accomplish the MC simulations with the MO-based PE surfaces in a practical central processing unit (CPU) time, the MO-based PE of each conformer was pre-calculated and stored before the MC simulations, and then only referred to during the MC simulations. The activation energies for transitions calculated by the MC simulations agreed well with the experimental ΔG determined by the NMR experiments. The analysis of the transition trajectories of the MC simulations revealed that the transition occurred not only through the transition states, but also through many different transition paths. Our computational methods gave us quantitative estimates of atropisomerism of the thienotriazolodiazepine compounds in a practical period of time, and the method could be applicable for other slow-dynamics phenomena that cannot be investigated by other atomistic simulations.
Stereoisomers often show different biological activity, pharmacokinetics and toxicity profiles because a ligand is three-dimensionally recognized by biological receptors. Therefore regulatory guidelines, such as Food and Drug Administration (FDA) guidelines,1) often request information on the biological properties of each stereoisomer in the drug development process. In addition to classical stereoisomers with chiral centers, much attention has been paid recently to the atropisomer.2–4) A compound which has axial chiralities caused by hindered rotations about rotatable bonds shows atropisomerism. Compounds showing atropisomerism have long-lived (1000 s or longer) conformers that can be isolated,5) and, like classical stereoisomers, each long-lived conformer often shows different biological activities.2,3)
Atropisomerism is caused by the high energy barrier between long-lived conformers due to the steric and/or electronic repulsions and is often observed in compounds with bulky substituents in the ortho positions. Half lives (t1/2) of the atropisomers are dependent on the energy barriers, thus the degree of atropisomerism is time-dependent.2–5) For compounds with an energy barrier of 120 kJ/mol or higher, the atropisomers have a t1/2 of several years at room temperature and behave like classical stereoisomers, but when the energy barrier is lower than 80 kJ/mol, the t1/2 is less than 10 s and compounds behave as mixtures or are simply flexible compounds. Compounds with a moderate energy barrier of ca. 100 kJ/mol have several hours of t1/2 at room temperature. From the viewpoint of quality control of pharmaceutical products and evaluating their biological properties, compounds with a moderate energy barrier are troublesome because racemization can occur within several hours at room temperature.
Estimating the energy barrier between atropisomers in the early stage of drug discovery and development is preferred because it will significantly affect the overall process of drug development. To evaluate the energy barriers around one acyclic rotational bond, LaPlante et al. applied a molecular orbital (MO) method and obtained satisfactory results.4) In addition to an atropisomerism caused by the rotational barrier of an acyclic bond, another kind of atropisomerism exists in the case of compounds with a medium-sized flexible ring. For such compounds, flexibility of the ring (ring inversion) causes atropisomers. Calculating the energy barrier along the transition paths for one acyclic chemical bond, in the manner of LaPlante et al., is straightforward and can be done by changing the torsion angle of the acyclic bond systematically and searching for the transition state (TS) of the compound. However, as LaPlante et al. pointed out, it is difficult to apply their method to estimate the energy barrier of compounds with a medium-sized flexible ring, where several torsion angles of rotatable bonds vary simultaneously and in concert.
Gilman et al. reported estimating the energy barriers of benzodiazepine analogs, which have a 7-membered ring.6) In their study, two torsion angles of the 7-membered ring of benzodiazepine were considered to define a single transition path between the atrpisomers and the energy barrier along the single transition path was calculated using the MM2 force field potential. In other studies not limited to compounds with a medium-sized flexible ring, the energy barrier between atropisomers was estimated by calculating the energy difference between the lowest energy conformer and the transition conformer by assuming a single transition path or a single transition conformer.7–10) In order to estimate the energy barriers between atropisomers with a medium-sized ring, which have some freedom of motion, it is preferable to consider all of the flexible torsion angles and the multiple transition paths across the energy barrier between the atropisomers, because the shape of the high dimensional potential energy (PE) surface for ring inversion is complicated and there may be more than one transition path to get across the energy barrier.
In this study, a computational method to estimate the energy barriers of transitions between atropisomers with a medium-sized flexible ring was developed and was applied to two thienotriazolodiazepine (TTDZ) analogs (Fig. 1), which have a 7-membered ring (all of the flexible torsion angles and the multiple transition paths were considered). In case that a single transition path of one acyclic rotatable bond is examined, a transition state of the compound can be identified straightforwardly as LaPlante reported. But, in case of multiple transition paths of the high dimensional PE surfaces with complicated shape such as a PE surface of medium-sized flexible, applying computational methods, such as Monte Carlo (MC) simulations or molecular dynamics (MD) simulations is required to seek many possible paths or trajectories in the high dimensional PE surfaces. In this study, the MC simulation method was selected so as to seek multiple transition paths, and the activation energies for transitions were directly estimated from the temperature dependency of the number of transitions observed during the MC simulations based on the Arrhenius equation. In order to make an MC simulation with a large number of MC steps (e.g. 1011 steps) possible, the PE surface, which controls transition dynamics between atropisomers, was calculated before the MC simulations, not during it, and the energy value of the pre-calculated PE surface for each conformer was simply referred to at each step of the MC simulation. A grid sampling of the whole PE surface using four torsion angles was performed and the potential energy for each grid point, which corresponds to each conformer, was evaluated at the HF/6-31G(d) level MO calculation. In order to validate the results obtained from our computational method, NMR experiments for the two TTDZ analogs were performed and the experimental energy barriers between the atropisomers were determined based on the NMR data and compared with the calculated energy barriers estimated by the MC simulations.
In this study, the atropisomerism of two TTDZ analogs, as shown Fig. 1, was evaluated computationally and experimentally as described below. These compounds possess a flexible 7-membered ring, thus they were considered to have a possibility of atropisomerism caused by the ring inversion of the 7-membered ring like other compounds with 7-membered rings.6,7,10–12)
Potential Energy SurfacesAs shown in Fig. 1, four rotatable torsion angles (ϕ1, ϕ2, ϕ3, and ϕ4) that cause the ring inversion of the 7-membered ring of TTDZ should be considered when evaluating the atropisomerism. As shown in Fig. 2, the lowest energy conformers and the PE surface for each compound were obtained by a grid search using the four torsion angles as follows. In the first step of the grid search for evaluating the PE surfaces, the initial conformers were generated in such a way as to rotate all of the four torsion angles systematically by 15° increments. The range of the four torsion angles was set as follows; 9 initial torsion angles ranging from −60 to 60° were sampled for ϕ1, 11 initial torsion angles from −75 to 75° for ϕ2 and ϕ3, and 24 initial torsion angles from 0 to 345° for ϕ4. Other conformers with ϕ1, ϕ2, ϕ3, and ϕ4 values outside the range listed above were ignored because such conformers form highly twisted 7-membered ring geometry and might not be passed during transitions between atropisomers. Based on the range and the increments defined, 26136 (9×11×11×24) initial conformers for each compound were generated and firstly optimized at the AM1 level by MOPAC200213) with the four torsion angles (ϕ1, ϕ2, ϕ3, and ϕ4) fixed. Among the 26136 optimizations, about 13000 energy optimized conformers with a valid 7-membered ring structure for each compound were obtained. Among the remaining optimizations, MOPAC2002 either gave no solution or gave conformers with a chemical structure different from the original structure. To calculate ΔE (which is the potential energy relative to the lowest energy for each compound) we needed to find the lowest energy for each compound; therefore, 162 (3×3×3×6) initial conformers were generated in such a way as to rotate all four torsion angles by increments of 60, and the lowest energy conformer at the AM1 level for each compound was obtained by the full optimization (no torsion angles fixed) of 162 conformers. From among the 13000 conformers generated for the PE surface of each compound at the AM1 level, we selected about 10000 conformers with ΔE of less than 80 kJ/mol, and these conformers were further optimized at the HF/6-31G(d) level by Gaussian0914) with the four torsion angles fixed. To obtain the lowest energy for each compound at the HF/6-31G(d) level, all of the low energy conformers that had been obtained in the full optimizations at the AM1 level, were fully optimized (no torsion angles fixed) at the HF/6-31G(d) level. The ΔE values at the HF/6-31G(d) level of about 10000 conformers for each compound were calculated and used to generate the PE surface, which is defined by the four torsion angles (ϕ1, ϕ2, ϕ3, and ϕ4) and ΔE. The ΔE value of the PE surface of about 3000 conformers, for which no optimizations at the HF/6-31G(d) level had been performed and no ΔE value was calculated, were simply assigned as a constant large value of 1000 kJ/mol, because these conformers might have an extremely high potential energy and are thought unlikely to participate in the transition processes. To summarize the process, the PE surfaces for each compound were generated by combining the information of the four torsion angles with the ΔE of 10000 conformers at the HF/6-31G(d) level and the ΔE of 3000 assigned as 1000 kJ/mol. These pre-calculated PE surfaces for each compound with the ΔE values at the HF/6-31G(d) level were used for the MC simulations.
One MC step of the MC simulation consisted of two procedures, a selection of a new conformer and a judgment to accept or reject this new conformer based on the potential energy of the new conformer. In the process of selecting a new conformer, a new conformer was randomly chosen from the neighboring conformers of the current conformer in the conformation space defined by the four torsions (ϕ1, ϕ2, ϕ3, and ϕ4), as described below. In this conformation space, the distance (Δϕij) between the conformers i and j were defined by the following equation (Eq. 1).
 | (1) |
where ϕk(i) was the k-th torsion angle (k=1, 2, 3, 4) of the i-th conformer. The neighboring conformers were defined as the conformers i and j with a distance Δϕij less than a constant value of Δϕc. In this study, Δϕc was set at 21.2 degrees which means that the conformational change at one MC step is limited to the change of only one or two torsion angles of 15°. Under these conditions, each conformer had an average of 25 neighboring conformers, and about 20% of the conformers had 32 neighboring conformers, which is the maximum number of neighboring conformers for one conformer.
To accept or reject a new conformer at one MC step simply obeyed the Metropolis method.15) If the ΔE of a new conformer was less than or equal to that of the current conformer, the new conformer was accepted. If the ΔE of the new conformer was greater than that of the current conformer, the new conformer was accepted if e−ΔΔE/RT was greater than a uniform random number between 0 and 1, where the ΔΔE was the ΔE of the new conformer minus the ΔE of the current conformer (ΔΔE=ΔEnew−ΔEcurrent), R was a gas constant and T was the absolute temperature. To find the ΔE for new and current conformers, the program simply referred to the ΔE value stored in the pre-calculated PE surface. These MC procedures were done using in-house programs.
Assuming the number of transitions between atropisomers observed during the MC simulations with the constant MC steps is proportional to the rate constant k, the activation energy, Ea, which corresponds to the energy barrier between atropisomers, could be calculated from the Arrhenius equations (Eqs. 2, 3),
 | (2) |
 | (3) |
where k is the rate constant of the transitions between atropisomers, ntrans is the number of transitions between atropisomers during the MC simulations, R is a gas constant, T is an absolute temperature and A and A′ are arbitrary values. Therefore, the number of the transitions between atropisomers during the MC simulation was counted up. MC simulations for each compound at nine different temperatures ranging from 500 to 900 K in 50K increments were executed to estimate the temperature dependency of the number of transitions. At each of the nine temperatures, 10 MC simulations with different random seeds and different starting atropisomers were performed. One MC simulation consisted of 1011 MC steps, which was thought to be long enough to cover a slow transition between atropisomers.
NMR ExperimentsNMR spectra of compounds 1 and 2 were measured in dimethyl sulfoxide (DMSO) solution at different temperatures to determine the coalescence temperatures (Tc) and the chemical shift differences in a slow exchange (Δυ0) of the methylene protons at position 6 (Fig. 1). At the coalescence temperature, the rate constant (k) of exchange between the atropisomers is defined as an Eq. 4.5)
 | (4) |
The Gibbs free energy (ΔG) of exchange between the atropisomers is calculated by the Erying–Polanyi equation (Eq. 5),
 | (5) |
where R is a gas constant, h is a Planck’s constant and kb is a Boltzmann constant.
For each compound, the two lowest energy conformers, aS_min and aR_min, were identified after the sequence of full optimizations at the AM1 level starting from 162 initial conformers followed by the full optimizations at the HF/6-31G(d) level starting from low energy conformers at the AM1 level, according to the process shown in Fig. 2. Three dimensional (3D) structures of aS_min and aR_min for each compound are shown in Figs. 3a and b. The lowest energy conformer aS_min was a complete mirror image of aR_min in both compounds 1 and 2. The crystal structure of compound 1 determined in-house superimposed well with the structure of the lowest energy conformer aS_min of compound 1 at the HF/6-31G(d) level, as shown in Fig. 3c. The starting conformers of the MC simulation, aS_MC and aR_MC, were defined as the conformers with four fixed torsion angles (ϕ1, ϕ2, ϕ3, and ϕ4) that were closest to those of the lowest energy conformers aS_min and aR_min. The details of aS_MC and aR_MC are listed in Table 1. The number of the transitions between aS_MC and aR_MC observed during the MC simulations was counted up.
Carbon atoms are colored in green, nitrogen atoms in blue, sulfur atoms in yellow, chloride atoms in dark green and hydrogen atoms in white. c) Superimposition of the lowest energy conformer aS_min and the crystal structure of compound 1. Carbon atoms of the lowest energy conformer are colored in green, carbon atoms of the crystal structure are in magenta, other atom colors are the same as a) and b). Hydrogen atoms are not depicted. 3D coordinates of these lowest energy conformers are listed in the supplemental material.
| ϕ1 | ϕ2 | ϕ3 | ϕ4 | ΔE (kJ/mol) | PC-1 | PC-2 | ||
|---|---|---|---|---|---|---|---|---|
| Compound 1 | aS_min | 47 | −66 | 63 | −159 | 0 | ||
| aS_MC | 45 | −60 | 60 | −165 | 2.2 | 93 | 23 | |
| aR_min | −47 | 66 | −63 | 159 | 0 | |||
| aR_MC | −45 | 60 | −60 | 165 | 2.2 | −93 | −23 | |
| Compound 2 | aS_min | 42 | −65 | 63 | −162 | 0 | ||
| aS_MC | 45 | −60 | 60 | −165 | 1.2 | 93 | 24 | |
| aR_min | −42 | 65 | −63 | 162 | 0 | |||
| aR_MC | −45 | 60 | −60 | 165 | 1.2 | −93 | −24 |
The PE surfaces at the HF/6-31G(d) level for compounds 1 and 2 were calculated as described in the Experimental section. Rough images of the shape of these PE surfaces are shown in Fig. 4. To create these two-dimensional PE surface maps, the dimensionality of the PE surface, defined by the four torsion angle values (ϕ1, ϕ2, ϕ3, and ϕ4) and one potential energy value ΔE, was reduced through two reduction steps. Firstly, the information of the ϕ4 angle was reduced because ϕ4 is the torsion angle that defines the rotation of the thiophene ring attached at 4 position of TTDZ core and the conformation of the 7-membered ring of TTDZ core itself can only be defined by the torsion angles ϕ1, ϕ2, and ϕ3. Among the set of conformers that have the same combination of ϕ1, ϕ2, and ϕ3 angles, but have a different ϕ4 angle, the lowest energy conformer was selected as a representative conformer for each combination of ϕ1, ϕ2, and ϕ3. Secondly, the information of the three torsion angles ϕ1, ϕ2, and ϕ3 was reduced to two dimensions by using principal component analysis (PCA). The conformers selected in the previous step were projected onto the two-dimensional map using the first and second principal components (PCs) derived from the PCA of the three torsion angles ϕ1, ϕ2, and ϕ3 (Fig. 4). In the two-dimensional PE surface map, the ΔE was represented by the color. The contribution ratios of the first and second PCs were 0.59 and 0.32 for compound 1 and 0.60 and 0.32 for compound 2, respectively. The high contribution ratios of the first two PCs suggest that there mainly exist only one or two conformational freedoms for the 3D structure of 7-membered rings of compounds 1 and 2 and they probably stem from the geometrical restrictions of 7-membered ring formation16,17) and three planar non-rotatable bonds (one double bond and two aromatic bonds) in 7-membered rings.
Each square denotes a representative conformer of each combination of ϕ1, ϕ2, and ϕ3 (see the main text for more detail). The ΔE from 0 to 20 kJ/mol is colored in cyan, 20–40 in blue, 40–60 in green, 60–80 in brown, 80–100 in magenta and >100 kJ/mol in red.
The existence of many red points, which mean high PE conformers, in Fig. 4a clearly show that the energy barriers between aS and aR for compound 1 were much higher than those for compound 2. In addition to the number of high energy conformers, Fig. 4b suggests that a transition between aS and aR without visiting any red-colored points might be possible in the case of compound 2, but in the case of compound 1, the transition without passing through high PE conformers in red might be difficult because most of the regions between aS and aR were covered by red points and aS and aR were separated by the conformers colored in red (Fig. 4a). The high energy barriers of compound 1 shown in Fig. 4a were caused by the steric hindrance between the methyl group at position 3 of the TTDZ core and the thiophene ring at 4 position. The steric effect of the thiophene ring at 4 position was supported by the results of our MC simulations which showed that in the 4H-TTDZ compound, which is a TTDZ analog with the hydrogen at 4 position instead of the thiophene group, the difference of the estimated Ea with methyl group at 3 position and without the methyl group becomes much smaller compared with that between compounds 1 and 2 (data not shown).
In order to confirm that the selection of the conformers in the first step by the AM1 level calculations for further calculation of the PE surfaces at the HF/6-31G(d) level (Fig. 2) is valid, we compared the ΔEs (obtained when generating the PE surfaces) between the conformers with four identical torsion angles at two different levels. The equations, ΔE(HF/6-31G(d))=1.96 ΔE(AM1) (r2=0.75) for compound 1 and ΔE(HF/6-31G(d))=1.84 ΔE(AM1) (r2=0.64) for compound 2, were obtained from simple linear regression analyses of the ΔEs at both levels. These regression models suggested that the ΔE at the HF/6-31G(d) and at the AM1 levels were linearly correlated and that the selection of the conformers by the AM1 level calculations in the first step was valid. These regression models also suggested that the absolute value of the ΔE at the HF/6-31G(d) level was about 2 times greater than the ΔE at the AM1 level, meaning that the PE surface at the AM1 level might underestimate the energy barriers. For this reason, there existed conformers with ΔE at the HF/6-31G(d) level greater than 80 kJ/mol (magenta or red squares) in the two-dimensional PE surface maps (Fig. 4), even though only conformers with ΔE less than 80 kJ/mol at the AM1 level had been selected for optimization at the HF/6-31G(d) level when generating the PE surfaces.
Monte Carlo SimulationsThe Arrhenius plots of the number of the transitions between aS_MC and aR_MC observed during the MC simulations with 1011 MC steps are shown in Fig. 5. The number of transitions observed during the ten MC simulations for each temperature was averaged, and the relation between the inverse of the temperature and the natural logarithm of the averaged number of the transitions was plotted. In the case of compound 1, the number of transitions at 500 K and 550 K were excluded from the analysis because the transitions occurred rarely (0–5 times at 500 K and 6–12 times at 550 K) during the MC simulations with 1011 steps and, from the statistical point of view, the number of transitions estimated from these rare observations was less reliable. The activation energy Ea for transition between atropisomers was estimated by the Eq. 3, ntrans=A′e−Ea/RT because a strong linearity between the natural logarithm of the number of transitions and the inverse of the absolute temperature (r2∼1) in both compounds 1 and 2 was observed, as shown in Fig. 5. The Ea values of compounds 1 and 2 were estimated to be 89 kJ/mol and 60 kJ/mol, respectively.
The number of transitions over the ten MC simulations at each temperature was averaged. For compound 1, the number of transitions at 500 K and 550 K are not shown (see the main text).
When the MC sampling details such as the Δϕc value were changed, the estimated Ea values were almost identical, and strong linearity (r2∼1) between the natural logarithm of the number of transitions during the MC simulations and the inverse of the absolute temperature (Arrhenius type relation) were constantly observed. For example, in the case of compound 2, the estimated Ea value was 61 kJ/mol when Δϕc was changed to 15, which meant only one torsion angle was allowed to change in one MC step, and the Ea was estimated to be 60 kJ/mol when Δϕc was changed to 26, which meant up to three torsion angles were allowed to change in one MC step, while the Ea was 60 kJ/mol when Δϕc was 21.2, as described in the Experimental section.
To identify the transition state (TS) conformers of both compounds, optimizations to find TS at the HF/6-31G(d) level were performed. The initial conformers for TS optimization were selected as follows: firstly, from among the trajectories at 600 K in ten runs of MC simulation, each transition trajectory between aS_MC and aR_MC was extracted. Secondly, for each transition trajectory, the highest energy conformer in the trajectory was selected. These highest energy conformers from the transition trajectories were collected and were used as initial conformers for TS optimization at the HF/6-31G(d) level. As shown in Fig. 6, two TS conformers for aS – aR transition with an almost planar 7-membered ring (ϕ1, ϕ2, and ϕ3 ∼0) were obtained for both compounds 1 and 2. It was confirmed that both TS conformers satisfied the necessary condition of a TS in having only one imaginary frequency that can be calculated by a normal mode analysis. The ΔE of these TS conformers were 94 kJ/mol for compound 1 and 65 kJ/mol for compound 2. These TS conformers are located at almost the center of the two-dimensional PE surface maps because the first and second PC values of these TS conformers were almost 0.
The TS conformers for each compound were mirror images of each other. The torsion angles (ϕ1, ϕ2, ϕ3, ϕ4) of TSs of compound 1 are (−2, −6, 9, 104) and (2, 6, −9, −104), and those of compound 2 are (−2, −9, 12, 123) and (2, 9, −12, −123). 3D coordinates of these TS conformers are listed in the supplemental material.
In order to examine whether the aS – aR transitions pass through the TS conformers or not, all of the trajectories of the transition between aS_MC and aR_MC observed during the MC simulations at 600 K and 800 K were traced. Transitions that did not pass through the TS conformers and their neighboring conformers (as defined in the Experimental section) were observed in about 40% of successful transitions of the MC simulation at 600 K and about 50% at 800 K of all of the transition trajectories for both compounds. This suggested that there exist paths for getting across the energy barriers between the atropisomers aS and aR without passing through the TS conformers.
To gain more information about the existence of the multiple transition paths, all of the transition trajectories observed in 10 MC simulations of compound 1 at 600 K were examined further. Among the 715 transition trajectories between aS_MC and aR_MC at 600 K, three representative trajectories are shown in Fig. 7. In case of the trajectory colored in red, compound 1 passed through the farthest point from the center, which corresponded to the TS conformers, of the two-dimensional map, and one other trajectory with a similar transition pattern was found from among the 715 transition trajectories. In the case of the trajectory in orange, compound 1 passed through the second farthest point from the center of the two-dimensional map, and three other trajectories with a similar transition pattern were found. In the case of the blue trajectory, compound 1 passed through the center of the map, where the TS conformer is located, with the smallest number of MC steps. Our method using the grid-based PE surfaces suggests that there are multiple transition paths between atropisomers, including ones not passing through TS conformers, at least at high temperatures.
The required CPU time for calculating the ΔE at the HF/6-31G(d) level per conformer was about two hours on average using one core of a Xeon X5640 CPU (2.93 GHz). Therefore, a total of about 20000 h (=2 h×10000 conformers) per compound is needed to calculate the PE surface at the HF/6-31G(d) level if one core of Xeon CPU is used. About 1100 h (=20000 h/18 cores, ca. 46 d) would be needed to calculate the PE surface for one compound if 18 CPU cores were fully used in parallel. One of the advantages of our ΔE calculations is that the ΔE calculation for each conformer is independent of others and so can be easily parallelized for the available number of CPU cores, without any special implementation. The required CPU time could be halved because the PE surfaces of compounds 1 and 2 had two-fold symmetry, which means that all of the conformers except two conformers ((ϕ1, ϕ2, ϕ3, ϕ4)=(0, 0, 0, 0), (ϕ1, ϕ2, ϕ3, ϕ4)=(0, 0, 0, 180)) had mirror-imaged conformers and thus had identical ΔEs.
The MC simulation with 1011 steps took about 20 h using one core of a Xeon CPU. A total of 1800 h (=20 h×9 temperatures×10 MC simulations, ca. 75 d) would be required to estimate the Ea of one compound if only one core of Xeon CPU is used. About 100 h (=1800 h/18 cores, ca. 4 d) would be taken to calculate the Ea for one compound if 18 CPU cores were fully used in parallel. In the same way as the ΔE calculations, the MC simulations are also independent of each other and can be easily parallelized. During the MC simulation by our method, the CPU was mostly occupied with generating the random number because no ΔE calculation is required at each MC step, since the pre-calculated ΔE value of the PE surface was simply referred to. If the energy of the conformer is evaluated at each step (as is usually done in MC simulations) and if the energy for each step is evaluated by a single point MO calculation at the HF/6-31G(d) level and if a single point MO calculation takes 0.004 h, then it would take 4×108 h (=0.004 h×1011steps, 46000 years) for one MC simulation of 1011 steps. In fact, 0.07 h was required for the single point MO calculations at the HF/6-31G(d) level for the compounds in this study in one CPU core and, therefore, we expect that ideally 0.004 (=0.07/18) hours per single point MO calculation would be required if 18 CPU cores were fully used in parallel. The difference between the duration needed for the MC simulation by our method and the estimated time required for MC simulations with energy calculations at each MC step demonstrates the power of the using the pre-calculated PE surface.
NMR ExperimentsIn Fig. 8, the NMR spectra of the methylene protons at 6 position of the TTDZ core of compounds 1 and 2 are shown. For compound 2, the slow exchange state was observed at 288 K and the chemical shift difference in slow exchange Δυ0 was determined to be 526 Hz. As the temperature was raised, an intermediate state appeared at 333 K and a coalescence state was observed at 338 K (Fig. 8b) and the coalescence temperature Tc of compound 2 was determined to be 338 K. At 343 K, compound 2 reached the fast exchange state. Based on the Eq. 4 and 5 and the values of 526 Hz of Δυ0 and 338 K of Tc, the ΔG of exchange between the atropisomers for compound 2 was estimated to be 63 kJ/mol. For compound 1, the slow exchange state was observed at 295 K and the Δυ0 value was observed to be 410 Hz. At 423 K, compound 1 was in the intermediate state but could not reach the coalescence state. Further elevation of the temperature was impossible because of a limitation in the hardware of the NMR device. Thus, the ΔG of compound 1 was estimated to be higher than 81 kJ/mol, based on the values of 410 Hz of Δυ0 at 295 K and >423 K of Tc.
Comparison of Calculated Ea and Experimentally Determined ΔGAs written in many text books, the following relationship between Ea and ΔG exists. From the Arrhenius equation, rate constant (k) is written as the Eq. 2.
| (2) |
From the Erying–Polanyi equation, k is also written as
 | (6) |
where kb is Boltzmann constant, h is Planck’s constant, R is gas constant and T is absolute temperature. Taking a natural logarithm of these Eqs. 2 and 6
 | (7) |
Since ΔG=ΔH−TΔS
 | (8) |
Differentiation by T
 | (9) |
Then Ea is obtained as
 | (10) |
Assuming the TΔS and the RT are small compared with the ΔG and are negligible, a comparison between the Ea calculated by the MC simulations and the experimentally determined ΔG of exchange between the atropisomers is possible. For compound 2, the activation energy Ea for the transition between atropisomers calculated by the MC simulations was 60 kJ/mol, which was almost equal to the experimentally determined ΔG, 63 kJ/mol. For compound 1, the calculated Ea of 89 kJ/mol was consistent with the experimental ΔG estimated to be greater than 81 kJ/mol.
Half-lives (t1/2) of atropisomers can be estimated from an equation written below
 | (11) |
where k is the rate constant, which can be estimated by the Erying–Polanyi equation (Eq. 6). The half-life of atropisomer t1/2 of compound 2 at 300 K was 0.01 s from the experimentally determined ΔG and 0.003 s based on the calculated Ea. These calculated and experimental t1/2 values suggested that the exchange of compound 2 between the atropisomers aS and aR is very fast and compound 2 will not display atropisomerism if the criterion for atropisomerism of having a t1/2 of more than 1000 s is applied.5) For compound 1, t1/2 at 300 K was experimentally determined to be greater than 14 s and t1/2 at 300 K was computationally estimated to be 350 s. The computationally estimated t1/2 of 350 s was used to judge the atropisomerism of compound 1 because the experimentally determined value, >14 s, of t1/2 includes both ranges of less than and more than the criterion, <1000 s and >1000 s respectively. It was suggested that compound 1 with t1/2 of 350 s will not exhibit atropisomerism, although the value of 350 s of compound 1 is numerically very close to the criterion, 1000 s.
It was interesting that the values of the potential energy ΔE of the TS conformers of compound 1 (94 kJ/mol) and 2 (65 kJ/mol) obtained from the analysis of the transition state were relatively close to the values of the activation energy Ea of compound 1 (89 kJ/mol) and 2 (60 kJ/mol) obtained from the MC simulations respectively. Although an analysis of the MC transition trajectories suggested the existence (at least, at high temperatures) of multiple transition paths and of transition paths that do not pass through the TS conformers, the relatively small difference between the ΔE of the TS conformers and the Ea derived from the MC simulations suggested that, at low temperatures, the contributions to the activation energy Ea of the transition paths that were far from the path along the TS conformers might not be significant. In other words, based on the MC simulations on grid-based PE surfaces in this study, there is a possibility that, at low temperatures, the compounds in this study get across the energy barriers between the atropisomers mainly at conformers near the TS conformers, even though there exist multiple transition paths at high temperatures. To verify the hypothesis described above, MC simulations at low temperatures with larger steps than 1011 might be needed, because the number of transitions between the atropisomers dramatically decreased at low temperatures.
In this study, a simple computational method of calculating Ea of transitions between atropisomers of a flexible 7-membered ring was proposed and was applied to the two TTDZ analogs. In our method, the ΔE of all of the possible conformers related to the transition between atropisomers was calculated at the HF/6-31G(d) level before the Monte Carlo (MC) simulations. Instead of calculating the ΔE at each MC step, as is often done for MC simulations, ΔE of the conformer at each step of the MC simulation was simply referred to the value of the pre-calculated PE surface. This procedure enabled an MC simulation with MC steps large enough to be applied to a slow transition phenomenon, such as an atropisomerism of a 7-membered ring with a high energy barrier, and 1011 MC steps were assigned for one run of the MC simulation in our case. Based on the Arrhenius equation, the Ea required for the transition between the atropisomers was estimated by the temperature dependency of the number of the transitions between atropisomers observed during the MC simulations. The calculated Ea for compound 2 was 60 kJ/mol, and this value agreed well with the ΔG for the transition between atropisomers, 63 kJ/mol, as determined by an NMR experiment. For compound 1, the experimental ΔG was estimated to be a range of value, >81 kJ/mol, rather than a specific value, because of a limitation of the NMR experiment. The Ea of compound 1 calculated by the MC simulations was 89 kJ/mol and was within the range of the experimental ΔG. The half-lives (t1/2) of the exchange between the atropisomers for compounds 1 and 2 were estimated to be 0.003 s and 350 s, respectively, based on the Ea values calculated by the MC simulations. It was suggested that neither compound will exhibit atropisomerism, based on the criterion that a t1/2 longer than 1000 s indicates atropisomerism. The analysis of the transition trajectories between atropisomers observed during the MC simulations revealed that, especially at high temperatures, the transitions occurred through multiple paths and there existed transition paths that did not pass through the transition state conformers. About 46 d of computation time were taken to calculate the PE surface of one compound at the HF/6-31G(d) level, which corresponds to about 10000 optimizations, with 18 cores of Xeon CPU used in parallel. It took about 4 d for the MC simulations of one compound, which corresponds to the 90 MC simulations with 1011 steps, using 18 cores in parallel. These computation times were very much shorter than the 46000 years that are estimated as required for one run of MC simulation at the HF/6-31G(d) level of 1011 steps when the potential energy is calculated at each MC step. Our method enabled MC simulations with an ab initio molecular orbital method-based PE surface and with MC steps large enough for slow transition phenomena to be run within a realistic computation time and the results agreed well with the experimental results. Further applications of this method to other phenomena with slow dynamics are expected.
