|Home | About | Journals | Submit | Contact Us | Français|
Incorporation of enediynes into anticancer drugs remains an intriguing yet elusive strategy for the design of therapeutically active agents. Density functional theory was used to locate reactants, products, and transition states along the Bergman cyclization pathways connecting enediynes to reactive para-biradicals. Sum method correction to low-level calculations confirmed B3LYP/6–31G(d,p) as the method of choice in investigating enediynes. Herein described as MI:Sum, calculated reaction enthalpies differed from experiment by an average of 2.1 kcal·mol−1 (mean unsigned error). A combination of strain energy released across the reaction coordinate and the critical intramolecular distance between reacting diynes explains reactivity differences. Where experimental and calculated barrier heights are in disagreement, higher level multi-reference treatment of the enediynes confirms lower level estimates. Previous work concerning the chemically reactive fragment of esperamcin, MTC, is expanded to our model system MTC2.
Esperamicin A1 (Figure 1) is a member of a larger family of naturally occurring products that share a novel molecular architecture and potent biological activity. Other family members include: dynemicin A, neocarzinostatin and calicheamicin γ1. The esperamicins were discovered in 1985,1 isolated from a fermentation broth of Actinomadura verrucosospora, and two years later their structures were reported.2,3 The enediyne compounds have drawn significant interest as potential sources of anticancer therapeutics because of their ability to cleave DNA, presumably through a Bergman cyclization.4 These compounds are among the world’s most powerful anti-tumor agents, as they exhibit potent activity against a variety of murine tumor models at injected doses in the 0.1 µg·kg−1 range.5,6 However, the enediyne natural products are not selective in their activity, and will cleave DNA in both healthy and cancerous cells. This indiscriminant behavior has spurred a flurry of experimental7–29 and computational17,30–69 research towards rational drug design. As of early 2008, only esperatrucin (Bristol-Meyers-Squibb) was in Phase II clinical trials, and two additional enediynes can be categorized as pre-clinical. A handful of enediynes are undergoing biological testing at other institutions. Promise for this compound class is found in the marketed drug Mylotarg (Wyeth) which is composed of calicheamicin tethered to a monoclonal antibody targeting acute myeloid leukemia.70,71
The exact mechanism responsible for esperamicin A1’s antitumor activity is not thoroughly understood. Several possibilities are presented in a recent paper by Capitani, et al.72 In the most commonly accepted mechanism for esperamicin’s activity, a bioreductive cleavage of the trisulfide tail is carried out by a reducing agent generating a thiol.73 The thiol then undergoes an intramolecular addition to the 10-membered ring containing the enediyne moiety, reducing differential strain between the reactant and transition state and increasing flexibility in the enediyne ring,74–76 which significantly lower the activation barrier for Bergman cyclization (Figure 2). This cyclization reaction is an electronic rearrangement where the 1,5-diyne-3-ene moiety rearranges to form a 1,4-didehydrobenzene diradical. Abstraction of hydrogen atoms from the DNA backbone by the diradical leads to single- and double-stranded DNA cleavage and eventual apoptosis. NMR studies of esperamicin and calicheamicin bound to DNA helices, as well as a computational docking study of dynemicin to DNA, provide structural rationale for the abstraction of hydrogen atoms from the DNA backbone.62,77,78
Recently, Alabugin et al. proposed an alternate mechanism for the role of the radical abstraction step in cycloaromatization kinetics.68 They argue for the feasibility of an intermediate intramolecular hydrogen abstraction step for benzannulated enediynes with methoxy substitution that would occur after Bergman cyclization and before the hydrogen abstraction from DNA. This intramolecular abstraction would essentially ‘quench’ the p-benzyne intermediate through the transfer of one of the radical centers to the methoxy group, making deactivation pathways such as the retro-Bergman cyclization inaccessible, effectively increasing the lifespan of the DNA damaging species.
Initial work involving the synthesis and investigation of enediyne fragments supported a relation between the critical internuclear distance of the forming carbon-carbon bond (Rcd; see Figure 2) and the cyclization barrier.79,80 This distance is dependent upon the acetylenic substitution and can be influenced by ring size if the acetylenes are connected through a larger macrocycle. Figure 3 includes the structures of many enediyne molecules (1–14) derived from the parent compound 1. Nicolaou et al. found that enediyne cyclization was spontaneous at distances below 3.20 Å, while above 3.31 Å the enediyne was stable at ambient temperature.79 It was concluded that the intervening length of 3.20–3.31 Å was critical for reactivity. This led to the design of molecule 11, a diol-containing analog of the parent enediyne 1, that has an internuclear distance of 3.20 Å and a measured activation energy of 31.5 kcal·mol−1.80 Compounds designed with larger internuclear distances were unreactive lending further support that this distance was the determining factor for reactivity. Schreiner argued for extension of this critical reactivity range to include internuclear distances spanning 2.9–3.4 Å after calculating the activation energies of monocyclic enediynes of varying size.35 In their analysis of enediyne reactivity, Alabugin and Manoharan have shown that the critical distance of 3.20 Å is indicative of spontaneous cyclization since the attractive two-electron interaction of the in-plane π-π* orbitals outweighs the repulsive contribution of the π-π interaction (likened to the transition-state of the symmetry-forbidden thermal [2s + 2s] cycloaddition).57
As more species were characterized, a competing explanation for the lowering of the cyclization activation barrier took the form of strain release across the reaction coordinate. In systems of reduced size such as 1–3, ring strain was not considered a decisive factor since a ring bridge such as that in the parent natural products does not exist. With the addition of fused ring systems, such as 12, reducing the ring strain in moving from the reactant to product became an equally important factor in rationalizing lower activation barriers.74,81
Increased reactivity for designed enediynes with high strain was found for a series of photochemically activated molecules.52 In this study, the internuclear distance ranged from 4.94–3.89 Å. Contrary to the critical distance hypothesis, the percent conversion at the longest distance was 100%, while at the shortest was only 66%. Clearly, both strain release and the distance between cycloaromatizing carbons play roles in determining the reactivity of enediynes.36,74,81,82 In addition to the compounds described above, extensive work characterizing other members of this family has been reviewed by Nicolaou and Dai.7
Designing promising new enediynes requires that the cyclization barrier be surmountable at body temperature. Once cyclized, the magnitude of the singlet-triplet (ST) gap controls reactivity.83 If the S–T gap is small, hydrogen abstraction via the triplet is quick since the singlet state is not overly stabilized (S-T gaps of ~2–5 kcal·mol−1 are indicative of reasonably reactive species). With a large gap, the triplet is destabilized relative to the singlet due to stronger through-bond interactions in the singlet, and the rate of abstraction is reduced. Early work at changing the electronic nature of the enediyne to alter the reactivity of the compound saw incorporation of a nitrogen atom into the double bond at position 3 which forms pyridine upon H-abs by the biradical.37,43,47,84,85 The S–T gap for the neutral species is high relative to 1, but protonation lowers this gap into a more reactive window (compounds 16–17). Design of an enediyne that contains a protonatable functionality is particularly appealing since selectivity to cancerous cells is related to cellular pH.86 The pH of a normal cell (7.5) drops to between 6.5–7.0 when cancerous, and can be lowered even further under hyperglycemic conditions. An ideal compound would be neutral and possess a large S–T gap (i.e. lower toxicity) when in the presence of healthy cells, but upon protonation in tumor cells its reactivity would resemble that of the natural products. Kraka and Cremer have utilized this design strategy to propose several new enediyne fragments incorporating amides and amidines (such as 18 and 19), but calculate retro-Bergman barriers lower than needed for functional drugs.87 Control of enediyne reactivity through pH has also been discussed in the density functional study of Alabugin and Monoharan.57 The competition between the retro-Bergmann reactions and H-abs are problematic for 16 and 17, as Cramer calculates retro-barriers to be inactivating.84
In their studies to investigate new enediynes with better drug-like properties, Kraka and Cremer have identified several qualifications that a designed molecule must meet, which include: 1) activity at body temperature, 2) ability of biradical to cleave DNA, 3) kinetic stability of biradical, and 4) low toxicity.43 These guidelines lead to molecules that upon cyclization do not decompose from competing reactions and have enthalpic barriers to cyclization less than 24 kcal·mol−1 at 310 K. Enediynes incorporated into rings commonly have barriers in this range. Initial work at computationally identification of novel enediyne replacements is promising. Even with this large body of work, considerable effort is still required before selective enediynes are identified. The problem of enediyne toxicity additionally remains challenging for the design of new drugs.71
Here we expand an effective model compound and evaluate methodology for the computational study of large enediyne fragments. Because of esperamicin A1‘s large size, and the rigorous theoretical treatment that its multi-reference electronic state demands, we have chosen to work with a series of smaller models that contain esperamicin A1’s chemically active core. The new model compound presented in this paper, MTC2 (Figure 3, structure 15), is a variation on the smaller MTC chemical model previously used by our group to assist in the study of the Bergman cyclization.88 The MTC2 model achieves a significant reduction in size by omitting the oligosaccharide arms that are important for esperamicin A1’s binding to DNA. Since we are not currently concerned with the triggering mechanism, the carbamate tail adjacent to the five-membered sulfur heterocycle has been pruned to further reduce the system size. These deletions reduce the calculation system size from one hundred sixty atoms to forty atoms, yet the MTC2 model incorporates esperamicin A1’s multi-cyclic ring structure, partially accounting for neighboring steric and electronic effects.
The Bergman cyclization is a challenging reaction to model because of the multi-configurational nature of the diradical species. Reliable description of the diradical state requires incorporation of both dynamic and static electron correlation effects. Several authors have published authoritative commentaries on the state-of-the-art in modeling these systems computationally.43,45,49,63,84,89 Multiconfigurational methods such as CCSD, CCSD(T), and CASSCF perform well for the parent enediyne system (e.g. compound 1), but are overly demanding for larger systems. Kraka and Cremer have shown CCSD(T)/pVDZ calculations perform well, and along with Grafenstein et al. have identified unrestricted broken spin DFT to account for the electron correlation effects needed to accurately model the reaction thermochemistry.45 Historical work performed with pure density functionals (e.g. BLYP) within the restricted formalism failed for modeling the enediyne cyclization (see discussion in Schreiner et al.63). Both the UHF-CC and RB-CC levels of theory have been shown to provide accurate diradical state descriptions of para-benzyne,54 but these methods are prohibitively expensive for all but the smallest enediyne model systems. Semiempirical theory has shown some promise, yet it is difficult to benchmark these methods without experimental results on the natural products themselves.33
Based on a wealth of theoretical calculations, Kraka and Cremer have identified B3LYP/6–31G(d,p) to perform adequately for geometries, followed by single-point calculations using larger basis sets, such as 6–311+G(3df,3pd), for providing accurate energies.45 Thus, in the current study, geometry optimizations of reactants (R) and product (P) triplets were performed using B3LYP90–93 with the 6–31G(d,p)94,95 basis set. Because of the open shell nature of both the diradical product (open-shell singlet) and, to a lesser extent, the transition-state species, all calculations on these structures were performed using BS-UB3LYP (broken-spin-symmetry, unrestricted) calculations. Both the “Guess=Mix” and “Stable=Opt” keywords were used to generate open shell singlet wavefunctions. When spin contamination was noted in the transition state, the transition-state triplet (single-point calculation for sum correction, see below) was computed.
A focused methodology was chosen for structures 1–30 consisting of low-level geometry optimizations followed by high-level energy evaluations, where the low and high levels are labeled MI (B3LYP/6–31G(d,p)//B3LYP6–31G(d,p)) and MII (B3LYP/6–311+G(3df,3pd)//B3LYP/6–31G(d,p)). For the new model compound (MTC2: 15), we obtained two sets of optimized geometries by using the B3LYP and the PBE1PBE96 density functionals with the 6–31G(d,p) basis set. Subsequent single-point calculations were computed with both DFT theories using the 6–311G(2d,p), 6–311+G(3pf,3dp), cc-pVDZ, aug-cc-pVDZ, cc-pVTZ and aug-cc-pVTZ basis sets.95,97–101 Additional single-point calculations on select structures used the CCSD(T)102–106 and BD(T)107 methodologies in conjunction with the aug-cc-pVDZ basis to address the accuracy of B3LYP energetics. Tight SCF convergence criteria were requested for all single-point calculations (10−8 Hartrees, or SCF=(Converg=8)).
Frequency calculations on all optimized structures allowed for determination of enthalpies and free energies at one atmosphere and at various temperatures. The reactant and diradical structures represent minima on the potential energy surface, and have no imaginary frequencies. The optimized transition-state structures have one imaginary frequency along the mode corresponding to bond formation in the Bergman cyclization. Computed thermodynamic properties include the reaction enthalpy, activation enthalpy, and activation free energy labeled ΔHr, ΔH‡ and ΔG‡, respectively. The retro-Bergman barrier is the difference of ΔH‡-ΔHr.
Triplet contamination of the product and transition-state singlets was removed using the sum method (both MI and MII) where the true singlet energy is backed out of the unrestricted singlet and triplet energy as weighted by the calculated spin contamination values (S2).45,101 The expectation value of the total spin operator is S2, where S2=s(s+1) and s is the quantum number for spin. Equation 1–Equation 4 describe the sum correction:
where E(UDFT,PS) is the energy of the product singlet using unrestricted formalism, E(PS) is the energy of the corrected product singlet being solved for, E(PT) is the product triplet energy (single-point calculation) at the contaminated singlet geometry, and S2 is the spin contamination of the Kohn-Sham orbitals (KS). The triplet energy is calculated at the geometry of the singlet, since contamination of the higher spin state is relevant at that specific geometry. The use of a sum formula was shown to dramatically improve results computed using B3LYP, but not for results calculated with pure functionals.43,45
Strain energies were calculated by first identifying portions of each molecule, or motifs, that remained constant (i.e., did not undergo a chemical change) throughout the reaction coordinate, and then truncating the molecule to that motif. Truncation of the R, TS and P geometries left open valences that were capped with hydrogen atoms. Minimizations were then performed on the two capping hydrogen atoms, keeping the rest of the motif constrained. Global minimizations were also performed on the capped motifs. The strain energy was then computed as Strain Energy = E(I of Scheme 1) − E(II of Scheme 1). This procedure is illustrated in Scheme 1, using reactant 2 as an example whose constant structural motif from R to TS to P is the phenyl ring. Note that atoms within the dotted areas indicate the minimized portion of the motif.
Absolute strain in reference to the global minima at each reaction coordinate, and changes in strain going from R to TS and TS to P were calculated. Strain was also estimated by performing MM2 minimizations of the reactants and products.108
Biradical stabilization energies (BSE) representing the enhanced stability of the biradical in comparison to the two possible mono-radicals was calculated according to Scheme 2.
Force field estimates of the strain energy were performed using the MM2 routine within Chem3D Ultra. All calculations were performed using revisions B.02 and C.02 of Gaussian 03.109
Structures of each system modeled (1–30) can be found in Figure 3–Figure 4. Experimentally determined reaction parameters for compounds theoretically characterized in this study can be found in Table I. The theoretical work concerning enediyne cyclization has been dominated by the characterization of 1 and a limited overview of the performance of methodologies is presented in Table II, including experimental values79,110,111. Pure density functionals perform poorly for this reaction while the hybrid functional B3LYP performs favorably. The best agreement between theory and experiment for 1 is achieved by performing single-point energy evaluation using extended basis sets, specifically 6–311+G(3df,3pd).
Biradical character in the product singlets and transition states is described in Table III. Spin contamination values spanned 0.94–0.99 for biradical P singlets, were 2.01 for every P triplet calculated (data not included), and ranged from 0.00 to 0.17 for transition states. The P singlet values are in the expected range of a 50:50 mixed representation of the singlet and no contamination values fell outside of the standard 10% variation indicating reliable energies. Spin contamination for a closed shell singlet transition state should be 0.0, and this was the common finding for the parent enediyne 1. When modeling larger systems, it was found that the contamination of the KS orbitals was substantial for several TS species as indicated in Table III, and indicative of the developing biradical nature for certain transition states. For this reason, the sum method was utilized to correct not only the P singlets but the relevant TS singlets as well. Correction factors found in Table III indicate that the sum correction to the P singlet (Sum:PS) lowers the energy of the singlet P. This correction factor is logically applied to the S–T gap, the retro-Bergman cyclization and toΔHr. The forward reaction barriers are affected by the sum correction to the TS singlet (Sum:TS,S) where the barriers are lowered by the indicated correction factors in Table III.
Table IV reports the M1, MII and MI:Sum and MII:Sum (sum corrected) enthalpy and free energy reaction barriers for each system that has known experimental values. The experimental values were determined at various temperatures, and the appropriate temperatures were used for each calculation. Summary statistics for each of the four methodologies are provided in Table V as mean signed error (MSE) and mean unsigned error (MUE).
Dipole moments of the reactants and transition states are provided in Table VI. As will be discussed in greater detail, two competing explanations for the lowering of the Bergman cyclization barrier are represented by the interatomic distance of the acetylenic termini (Rcd) and the molecular strain energy. Several measures of strain as well as Rcd can be found in Table VII and the relevance of each of these molecular descriptors are provided in Table VIII.
While a consistent methodology was applied to the investigation of structures 1–14 and 21–30, the novelty of structure 15 (MTC2) drove us to explore several complementary computational approaches, whose results are provided in Table IX.
Our analysis leads to the identification of MI:Sum as the most statistically accurate methodology for looking at enediyne cyclization for a diverse set of reactants. This is inconsistent with the literature findings for structure 1 alone. Table X provides a summary of the activation free energies of each species at 310 K calculated using MI:Sum. The CCSD and BD calculations used to confirm the lover level results on structures 6 and 7 are provided in Table XI.
Experimentally and theoretically determined values of ΔHr, ΔH‡, ΔG‡ and the S–T gap for 1 are found in Table I, ,IIII and III. As shown in Table II, B3LYP reaction enthalpies improve with relation to experiment as the basis set expands from double-zeta to triple-zeta (i.e. VDZ vs. TDZ or 6–31G vs. 6–311G). Reaction enthalpies are consistently underestimated for double-zeta calculations for 1; conversely, enthalpies calculated at the triple-zeta level are slightly overestimated. The B3LYP activation enthalpies show reverse trends with the double-zeta calculations predicting ΔH‡ accurately and the triple-zeta calculations overestimating the barrier by 4–6 kcal·mol−1. Both pure functionals consistently underestimate the barrier height, though expansion of the basis set from double- to triple-zeta reduces the errors from ~4 to ~1 kcal·mol−1. For a given methodology, accuracy in one experimental parameter is often associated with inaccuracy in another parameter, and choosing a methodology based off the results for one enediyne may be misleading (as has historically been done with 1).
Comparing the MI and MII calculations, MII consistently overestimates the reaction barriers, but closely predicts the reaction enthalpy. The S–T gap is underestimated by MI and MII by approximately 1.3 kcal·mol−1, while MI:Sum and MII:Sum overestimate the gap by approximately 1.4 kcal·mol−1. MI and MI:Sum best reproduce the experimental reaction barriers. Unfortunately, the MI:Sum reaction enthalpy is significantly underestimated by 7.8 kcal·mol−1. Basis set expansion leads to a better characterization of the biradical product but at the expense of computational cost.
Through an extensive exploration of computational methodologies, Kraka, Cremer and co-workers identified B3LYP/6–311+G(3df,3pd) single-point calculations on 6–31G(d,p) geometries as providing the lowest mean absolute error for electronic activation energies of the forward and reverse reactions and the electronic reaction energy for 1 (3 parameters).45 Furthermore, they found sum method corrected values were even more accurate with a MUE of only 1.2 kcal·mol−1. Their analysis relied on the comparison of electronic energies, and the correction of the experimentally derived enthalpies to electronic energies using computed vibrational correction factors. In this work we chose to compare the experimental ΔHr and ΔH‡ values directly with the theoretical values, which provided an MUE of between 1.9 and 3.8 kcal·mol−1 for each of the methods listed in Table II. Method MII has one of the lower error margins, in agreement with the analysis of Kraka, Cremer and co-workers.45 Table V’s error statistics for structure 1, using ΔHr, ΔH‡, ΔG‡ and the S–T gap parameters, indicate that MI and MII errors are nearly identical (MUE of 2.8 and 2.5 kcal·mol−1 for the sum corrected values). We find a larger associated error for MI and MII compared to the results of Kraka and Cremer for four reasons. First, we choose to include the S–T gap as a measurable experimental value for error diagnostics. Second, we choose to interpret the sum method corrections in a different way. Sum correction of the singlet product energy will alter the S–T gap, and the reaction enthalpy, but has no direct effect on the computed forward barrier heights. Kraka and Cremer rationalize that a more exothermic reaction (sum method lowering the energy of the product) will, by the Hammond postulate, lead to a transition state with a lower energy. Instead, we choose to use the product sum correction to correct the product energy only, and not the transition state energy; we altered the TS energy only if the TS showed biradical character as judged from the spin contamination expectation values. Third, we compared enthalpies and not electronic energies. Fourth, we chose not to include the retro-Bergman reaction as this is a mathematical combination of the reaction enthalpy and activation enthalpy that are already contained in the error function. That being said, for enediyne 1, we are in agreement with Kraka and Cremer that sum method correction improves the error statistics for this compound, and that MII provides an accurate treatment of the Bergman cyclization. Further computational support for the selection of MI and MII can be found in the work of Prall et al. where error diagnostics indicated B3LYP/6–31G* to be one of the methods in closest experimental agreement.89 Utilization of other functionals or more extended basis sets from the analysis of Prall et al. did not yield lower errors substantially enough level to warrant choosing an alternate methodology other than MI and MII used in this study.
The initial goal of this work was to investigate the reaction thermodynamics of a much larger enediyne system (15: MTC2) than previously completed. Since no experimentally determined parameters are available for this specific compound, it is necessary to ensure that our method is reliable. The computational literature cited in the introduction seldom strays from covering reactions other than the enediyne cyclization of structures 1–4. While lower level molecular mechanics, semi-empirical, or limited high-level calculations for other structures present in Table I exist, the selection of a suitable computational methodology for treating enediynes might depend upon the size and chemical diversity of the compounds being modeled. Given the extensive literature of experimentally determined parameters, we chose to further test methodologies MI and MII on compounds 2–15, and 21–30, which will be grouped together and discussed according to molecular similarity in the next sections.
Synthetic addition of a benzene group to the enediyne double bond was initially pursued to lower the cyclization activation barrier measured for 1 (ΔH‡ = 28.2 kcal·mol−1). Compound 2 is the most simplified version of a benzannulated species, and methoxylation of compound 2 provides 3. While the experimental S–T gap for 2 has not been determined,ΔHr,ΔH‡, and ΔG‡ are available. For 1, the spin contamination of the TS was zero resulting in the sum correction being applied only to the product energy. For 2, spin contamination is noted in the TS, as reported in Table III. The contamination is present only for MI and in moving to a more extended basis with MII this contamination disappears. The activation enthalpy calculated with MII is grossly overestimated although ΔHr and ΔG‡ are in better agreement with experiment than the corresponding values from MI (Table IV). Comparison can be made to the experimental results of Roth et al. as has been done in other studies, although other activation energy measurements of 2 in chlorobenzene indicate a consistent value for this parameter (an average value of 25.0±2.3 kcal·mol−1). For this reason, the Roth enthalpy was averaged with the other two measurements for the error analysis, indicating that MI performs better than MII.
The methoxylated version of 2, compound 3, was pursued as a proof-of-concept experiment to show the feasibility of intramolecular H-atom abstraction from an internal hydrogen source to the biradical intermediate.68 It was found that the methoxy group acted as a radical trap preventing the retro-Bergman, although the participation of such a mechanism in the natural products has not been established. While methoxylation alpha to the forming biradical serves to alter the secondary radical abstraction reactions, the initial reaction thermodynamics of cyclization are not altered given that the reported activation energies for 2 and 3 are within experimental error. The presence of the methoxy group lowers the S–T gap from 5.6 (2) to 5.0 kcal·mol−1 (3).
Compound 3 has two dominant conformations where the methoxy is pointed away from (anti) or towards (syn) the acetylene arms. At 310 K, the anti conformation is more stable than syn by a free energy difference of 4.4 kcal·mol−1 (MI), indicating that it is the dominant reactive conformation at this temperature based on a Boltzmann population analysis. Thus, all tabulated values are for the anti conformation even though the reaction barrier of the syn variant is calculated to be lower by 2.8 (MI) and 2.9 (MII) kcal·mol−1.
Zeidan et al. calculated the cyclization activation energy of 3 to be 24.7, 28.8, and 31.4 kcal·mol−1 with BS-UBLYP/6–31G*, BD(T)/cc-pVDZ//UBLYP, and BS-UB3LYP/6–31G*, respectively. While the pure functional was in better agreement with the reported experimental average of 25.2±1.5 kcal·mol−1 (Ea), the performance of pure functionals are often associated with higher error for other systems such as 1.45 The hybrid results are more in-line with those reported in Table IV indicating an overestimation of the computed reaction barrier (3) by 7.4 and 9.0 kcal·mol−1 for MI and MII. Sum method correction had no effect on these discrepancies as the TS species of 3 showed no biradical character. Neither MI nor MII is accurate for this species, but MI has less associated error.
Additional characterization of functionalized benzannulated compounds was performed by Zeidan et al., who examined the effect of aldehyde and nitro substituents.26 The authors found that steric compression (leading to a lowering of Rcd) of the enediyne by ortho substituents led to lowering of the calculated cyclization barriers.26 Compounds 23–25 were experimentally measured to have enthalpic barriers of 21.7, 20.7 and 23.6 kcal·mol−1. In all three cases, the barriers are lowered in comparison to the parent compound 2. This lowering was due to the ortho-effect described and not due to ring-closure. According to the Rcd hypothesis, shorter internuclear distances are indicative of lower barriers; for this set of compounds, and including 2 and 3 as references, Rcd is 3.94, 4.02, 4.16, 4.20 and 4.21 Å for 24, 23, 3, 25 and 2, respectively. The experimentally measured enthalpic barriers span the range 20.7–25.0 kcal·mol−1, and all reactants follow the Rcd trend except for 25, which falls between 23 and 3. The strain release hypothesis can be equally invoked since the build up in strain from the reactant to transition state steadily increases from 2.6 to 3.3 kcal·mol−1 for this series. When calculating barrier heights with MI, MI-sum, MII and MII-sum the reactivity trend switches order to become 23, 24, 2, 25 and 3, where the four methods to calculate enthalpic barrier heights are on average overestimated by 6.9, 8.3, 6.0, 7.5 and 8.2 kcal·mol−1, respectively. The error for MI:Sum (all compounds; Table V) is 3.0 kcal·mol−1, such that the average error of 5.7 kcal·mol−1 for this series of five compounds is disappointing compared with the study as a whole. For this series of compounds Zeidan et al. saw better agreement between measured activation energies (26.2, 21.3, 24.2 and 22.3 kcal·mol−1 for 2, 24, 25 and 23, respectively, Table I) and either BLYP (24.5, 20.9, 24.1 and 20.6 kcal·mol−1) or MP2 (25.3, 22.1, 24.9 and 22.1 kcal·mol−1) versus B3LYP (31.3, 27.9, 30.9 and 27.7 kcal·mol−1).26,68 Zeidan calculations utilized the 6–31G** basis set and do not appear to have been sum corrected. While the pure functional was not applied globally in the present study, it may be of use to return to the structures presented here without hybrid functionals. For the compounds where a benzene group has been annulated to the enediyne double bond (2, 3 and 23–25) there is no clear choice to which method is most suitable for thermodynamic characterization. However, Jones and Warner found BLYP to have a higher barrier than experimentally determined for 2, and attribute this discrepancy to its lowered retro-Bergman barrier in comparison with 1.69 Similar to the present study, compounds 1–2 cyclize with nearly identical calculated barriers suggesting that annulation in this case has little effect on the forward barrier height. Using the arguments of Jones and Warner, it may be that experimental results showing variation in the activation barriers because of benzannulation actually result from the retro-barriers.
Nicolaou et al. pursued the ring closure (macrocyclize) of 1 in an attempt to lower the cyclization activation barrier and arrived at compounds such as 4, 10 and 11. The experimental activation free energy for 4, 10 and 11 are 24.6, 26.8 and 24.8 kcal·mol−1, respectively, indicating tight agreement in the reaction free energy barriers (because the experimental temperature range was 333–378 K for 10, we chose the median temperature of 350 K for computing its thermodynamics). Computationally, the MI-sum method best reproduces the activation free energy, with values ranging from 24.4–25.2 kcal·mol−1 for 4, 10, and 11. Based off similar scaffolds where the pendant hydroxyl groups are distal to the reaction center, the similarity in reaction barriers (ΔG‡) is not surprising. However, there is less agreement for the reported enthalpic barriers.
Compounds 4 and 11 measured in benzene solution provide enthalpic barriers of 23.2 and 23.0 kcal·mol−1. However, a second measurement of 11 in THF provided a value of 30.9 kcal·mol−1, a surprising result considering the magnitude of the solvent’s dielectric constants are not significantly different. A possible explanation may lie in the π-stacking ability of the solvent since this property is vastly different in comparing THF to benzene. As proposed by Russell, an aromatic solvent such as benzene might stabilize the transition state through π-stacking of the developing aromatic ring (KC Russell, personal communication). Thus, reactions run in aromatic solvents might have lower activation parameters compared to those run in non-aromatic solvents or in comparison to gas-phase calculated values. Compound 10 was cyclized in THF, and it is not surprising that there is agreement in the reported enthalpy of 30.2 kcal·mol−1 and that of 30.9 kcal·mol−1 reported for 11. While the aromatic solvent hypothesis is interesting, it does not hold for compound 2, since the enthalpic barriers measured in methanol and the aromatic solvent chlorobenzene are nearly identical. However, there may be a normalizing effect of the higher dielectric constant of methanol. The rate enhancement at lower dielectrics was not large in magnitude for another study where Feng et al. measured the half-life for an aza-enediyne cyclization to vary from 3.6 to 8.1 minutes over a dielectric range of 1.89–37.5.112
Compound 27 contains the benzannulation component as well as the ring closure component. Benzannulation led to lowered cyclization barriers compared to ring closure. The effects do not appear to be additive, as the experimentally measured enthalpic barrier (24.4 kcal·mol−1) is in close agreement with the other members of the benzannulated family. Compared to the benzannulated species, the ring closure compounds all have shorter Rcd distances with reduced strain build up for 4 and 10. Strain release for 11 and 27 is evident as opposed to strain build up. Given the similarity of molecular composition, and the tight range of calculated MI:Sum enthalpic barrier heights (Table IV) of 23.3–25.5 kcal·mol−1, the reported enthalpic barriers of 30.2 and 30.9 kcal·mol−1 for 10 and 11 seem high. A second reported measurement of 23.0 kcal·mol−1 for 11 agrees better with the calculated enthalpic barrier. In summary, for the ring closure compounds, MI:Sum best reproduces experimental enthalpic and free energy barriers when the experiment (enthalpic barrier determination) occurred in aromatic solvent.
Russell and coworkers have synthesized and characterized several benzannulated enediynes with heteroatom substitution ranging from a single nitrogen (5) to nucleic acid type rings (6).113–115 Several studies investigated the effect of solvent and tautomerization on enediyne reactivity. In the case of 9, the reactivity was found to be dependent upon solvent polarity where moderately polar solvents such as acetonitrile and methanol induced much longer half lives compared with non-polar solvents.115 The dielectric constant of THF is higher than that of benzene, however the measured half-life of 9 is 96 minutes in benzene and only 16 minutes in THF contrary to the trend for moderately polar solvents. In comparing solvents with dielectrics on the order of 30 to those on the order of 2–7 (a greater spread), it was apparent that increased polarity raised the activation barrier, in line with Feng et al.112 This was attributed to greater polarity in the reactants as compared to the transition states, or in other words, increased solvent polarity stabilized the reactant to a greater extent leading to an increased reaction barrier height. Calculations on four of the five heteroatom structures (5–7 and 9) show the reactants’ dipole moment to be more polar than the transition state (see Table VI), supporting the solvent polarity theory. Compound 8 has increased polarity in the TS going against this trend. Since the polarity of reactant 9 is greater than the corresponding TS, it is likely that the proposed solvent polarity mediation of the reaction rate is playing a factor in slowing the reaction in methanol and acetonitrile.115
Russell et al. have provided evidence that the concentration of 1,4-cyclohexadiene (1,4-CHD) does not alter the reaction kinetics.113,115 Since 1,4-CHD acts as a hydrogen atom source in the determination of some reaction thermodynamics, the concentration of this reagent may influence the rate of reactant loss/product formation. The concentration of 1,4-CHD was varied over the range 0–600 nM and in the case of 8 the half-life for the four measurements varied with an average and relatively small error of 113.8±6.7 minutes. They point out that one would expect this variation to be greater, and in fact Zeidan et al. did find a rather large effect on kinetic rates for 2 and other ortho-substituted variants depending on 1,4-CHD concentration.26 At low concentrations of hydrogen-atom donor, they found the rate to be highly dependent on the 1,4-CHD concentration, but the rates stabilized at concentrations above 0.4 M. At low concentrations of 1,4-CHD (<0.1 M) the variation in reported activation energies could be as much as 4.8 kcal·mol−1.26
These findings have significant impact on the experimental design used to measure the activation parameters for compounds 6–9. At low concentrations of hydrogen atom donor, competing reactions for the loss of reactant may dominate. Contrary to this, when concentrations of hydrogen atom donor such as 1,4-CHD are in excess, kinetic measurement should be applicable only to Bergman cyclizations. Once a given reactant has undergone Bergman cyclization the product biradical must hydrogen atom abstract according to a second rate constant which must be faster than the competing reverse rate constant for ring opening. Not only this, but the reactant may undergo competing side reactions such as polymerization; such side reactions would effectively remove reactant, resulting in determination of an artificially low activation barrier to cyclization (e.g., polymerization of benzyne with enediyne).17,26,116–119 From the work of Zeidan and co-workers, the 1,4-CHD concentration should not be set below 0.2–0.3 M, far above the concentrations used for 6–9 (1.2×10−7 M 1,4-CHD). Alabugin and co-workers find that the threshold 1,4-CHD concentration is unique to each enediyne studied.26
When comparing the measured activation energies for an ortho-NO2 substituted version of 2, the Ea varied from 19.1–23.3 kcal·mol−1 over the 1,4-CHD range of 0.03–0.58 M, as measured through the disappearance of the reactant (Ea’s calculated by monitoring product formation were within experimental error).26 Zeidan et al. provide a clear breakdown of the reaction kinetics in comparing the pseudo-first order rate constants, as is assumed at varying concentrations of 1,4-CHD in the literature, and the true rate of the first-order cyclization constants. In decomposing the first-order rate constant for the loss of reactant 25, the activation energy is calculated to be 25.4±0.2 kcal·mol−1 compared to 23.9±0.3 kcal·mol−1 reported for the pseudo-first order measurement.
The parent benzannulated species, 2, has been experimentally characterized in methanol and chlorobenzene solvents. For this species, little difference was seen in the reported activation enthalpies as derived from the activation energies (Table I). Incorporation of a nitrogen into the benzene ring, compound 5, lowers the activation barrier by 4.1 kcal·mol−1 with respect to the average ΔH‡ of 25.0±2.3 kcal·mol−1 for 2. The solvent dielectric of CDCl3 (4.9) lies just outside the range of solvents used for the determination of 2 (5.4–33.0), so solvent variations should not be influencing barrier height differences to a large extent. Calculated values in Table IV show closer agreement between 2 and 5, so it is not clear that nitrogen incorporation outside of the quintessential Bergmann cyclization atoms (i.e. 1) play a large role in lowering the barrier.
When comparing structures 1 and 16, Cramer computed the activation barriers to be 27.7 and 19.4 kcal·mol−1, respectively.84 These values were calculated using a composite enthalpy derived from BCCD and CCSD calculations using the cc-pVDZ basis set, and show a substantial barrier lowering upon nitrogen incorporation (in this case, nitrogen substitution was a direct component of the forming 6-membered ring). Cramer found an increased biradical stabilization energy (BSE) for the aza-product (-14.9 for 16 vs. −4.5 kcal·mol−1 for 1) through interaction of the nitrogen lone pair and the p-biradical. This increased BSE led to higher exothermicity of the reaction for 16 and according to the Hammond postulate lowered the activation enthalpy, as reflected in the calculated barrier heights. Using sum correction to the B3LYP BSE, MI:Sum provides stabilization energies of −10.7 and −3.7 kcal·mol−1 for 16 and 1. Compound 5 does not have a nitrogen p-orbital alpha to the radical center to provide such a large BSE. In this case, the BSE for 5 is only −3.9 kcal·mol−1 at MI:Sum, and the corresponding value for benzannulated 2 is - 4.2 kcal·mol−1. Less computational demanding B3LYP BSE values are qualitatively in line with the CCSD(T) values although quantitative agreement suffers. Pursuing BSE at CCSD(T) for the structures reported here is computationally troublesome. Not only are the BSE values the same for 2 and 5, the critical bond distance and strain parameters are equivalent, complicating any attempts at rationalizing experimentally observed variation in reactivity of these two molecules. These two species cyclize with almost equivalent barrier heights at MI:Sum, and the retro-barriers do not vary significantly.
Incorporation of additional heteroaromaticity into the benzannulated ring as with 8 led Russell to measure one of the fastest rates of cylization for an enediyne with an experimental activation energy of only 16.8±0.9 kcal·mol−1.113,114 This represents a cyclization rate less than that of the natural product calicheamicin (activation free energy of 19.3±0.2 kcal·mol−1, Table I). Converting the reported activation energy to a free energy at 310 K using our calculated ΔS provides a free energy barrier of 18.1 kcal·mol−1 for 8. The measurement of 8 was performed in a non-aromatic non-polar solvent, and the rates of 5 and 9 were measured concurrently. Nitrogen incorporation was invoked as the cause of the barrier lowering when comparing 2 to enediynes 5 and 8.113,114 However, the secondary benzannulation, in moving from a structure of type 5 to 9, was speculated as the factor responsible for increasing the barrier height a substantial 12.1 kcal·mol−1 as judged by ΔEa. Further work by Russell et al. indicated that the tautomeric state of 8 had effected the rate of cyclization.113 Compounds 6–8 were measured under the same conditions and showed an increased barrier for the keto-version 6. Calculation of half-lives indicated compound 8 to have the longest t1/2 of 115.5 minutes at 120°C followed by 17.5 and 4.6 minutes for 7 and 6, respectively. Compound 8 had the longest half-life, but smallest cyclization barrier of the three compared. This seeming discrepancy, in that the reactivity of the series is not controlled by the activation parameters, was explained via two proposals. First, the ground state energies of 6 and 7 were lower than the lactim version such that fewer reactant 8 molecules passed over the reaction barrier (the tautomeric preference of 7 over 8 is 22.1 kcal·mol−1 using MI, supporting this notion). Secondly, the variation in observed half-life was attributed to entropic effects not captured in the determining the activation energies. Entropic changes (Table X) indicate that entropy effects associated with 8 lie intermediate to those of 6 and 7. Correction of activation energies to free energies at 298 K does not change the order of reactivity for 6–8.
From an activation energy perspective, compounds 6–8 have increasing half-lives (less reactivity) in an opposite trend to the activation barriers. The free energy barriers calculated with MI for 6, 7, and 8 are 30.7, 30.4 and 31.6 kcal·mol−1 at 310 K. These barrier heights correctly predict 8 to be least reactive, with a calculated barrier difference of 8 versus 7 to be 1.2 kcal·mol−1, leading to 7 cyclizing 7-fold faster at 310 K. This variation is not from increased entropic penalties as the entropy change at 310 K for 7 and 8 is calculated to be lower for the lactim version. The 7-fold rate change is in near agreement with compound 7 reacting 6.6-fold faster at 120°C compared to 8 as judged by a ratio of half-lives. However, it must be considered that the ratio of 6:7 is not reproduced by comparing free energies and that the calculated free energies are all within methodological error of each other (herein assumed to be the MUE of MI:Sum = 2.8 kcal·mol−1).
Turning to the interatomic distance Rcd, this distance increases from 4.11 to 4.22 to 4.29 Å for 7, 8, and 6, indicating that the lowest calculated activation free energy corresponds to the compound with the shortest critical bond distance, however the reactivity trend in half-lives does not follow Rcd, nor does is strictly follow the trend in build up in molecular strain from R to TS as judged by ΔTS strain (Table VII). Here, the strain built up in 8 is greatest followed closely by 6 and then 7. The reported half-life determinations are in better agreement with the theoretical values compared to the measured activation energies.
Of the structures studied by Russell, the most extreme activation barrier was found for 9, where this compound deviated from the others measured by up to 17 kcal·mol−1. Critical bond distances and molecular strain show 9 to have a bond distance within the range of the other compounds; similarly the strain value of 4.9 kcal·mol−1 implies a small variation in activation barrier compared with the rest of this class of enediyne. Stemming from these similarities and the agreement in calculated enthalpies, it is likely that compounds 5–9 should all cyclize with barriers similar to 9. Compounds 5–9 herein considered to be a sub-set of the whole have an average MI-sum ΔH‡ and ΔG‡ of 29.1±0.6 and 31.1±0.5 kcal·mol−1 indicating little variation in predicted reactivity. While variation in reactivity may be seen at different temperatures sampled during cyclization experiments, for biologically active agents cyclizing at 310 K, Bergman cyclization of compounds 5–9 would not be expected to occur.
Because of the lack of agreement between the above described compounds and the reported experimental values, we chose to confirm our findings for 6 and 7 using computationally demanding CCSD(T) methodology and just BD(T) for 6. While we are not trying to comment on benchmarking code and timings, we must point out that this calculation level was incredibly time demanding, and prohibited its use for investigating all structures. Comparing electronic barrier heights for compounds 6 and 7, values of which are found in Table XI, it can be concluded that the B3LYP values are overestimated by approximately 3.0 kcal·mol−1 for 6 using either multi-reference method as a reference point. BD(T) calculations were not pursued for 7 or any other species, as we think the data presented in Table XI supports the notion that the barriers for 6 and 7, and for that matter, the class of benzannulated species itself are accurately treated by B3LYP. While the multi-reference treatment brings the values into closer agreement with the divergent experimental values, they are still off by 6.2 and 7.4 kcal·mol−1 for 6 and 7, respectively (if B3LYP thermal corrections are applied to the CCSD(T) values).
The underestimation in the experimentally measured cyclization barriers for this sub-set is likely due to monitoring of reactant loss only and not product formation. Competing side-reactions lowering the concentration of reactant, in conjunction with low 1,4-CHD concentrations, likely led to low measured barriers.
After it was determined that connecting the acetylene arms through a closed ring system affected the rate of reaction, by bringing the otherwise large barrier to cyclization (e.g. 1) into more accessible ranges (i.e. enthalpic barrier <24 kcal·mol−1), compounds 12–14 were synthesized in order to better understand the interplay of critical bond distance (Rcd) and strain energy. In order to reduce the computational demand of these systems, the t-BuMe2Si- group, shown in Figure 3, was truncated to a methyl group; this substitution is not expected to play a significant role in the reaction thermochemistry.
Experimental activation free energies were determined to be 26.3, 27.4 and 32.0 kcal·mol−1 for 12, 14 and 13, respectively, at various temperatures (Table I and Table IV). General trends will be considered since all three compounds did not have thermodynamic parameters determined at the same temperature. What is clear is that compound 14 cyclizes fastest and contains a fused five-membered ring with a bridgehead sp3 center. Kinetic measurements of ring closure for structures 12–14 indicate ring strain to outweigh acetylenic proximity. The distance Rcd for compounds 12 and 13 from crystallographic analysis was found to be 3.391 Å and 3.368 Å, which fall just outside Nicolaou’s originally proposed range of reactivity. Cyclization reaction rates indicate that compound 13 reacts 650 times more slowly than 12 at 124°C even though Rcd is shorter (albeit by only 0.02 Å given the assumption that the oxime of the solved crystal structure in place of the carbonyl of 13 does not effect Rcd). When the sp2 center bridgehead keto is converted to an sp3 center alcohol, the rate of reaction is 216 times faster for 14 than for 13, and only 33% of 12’s rate. Calculation of strain relief using molecular mechanics allowed Magnus et al. to rationalize the rates of reaction as stemming more from strain relief than the separation of the acetylenic carbons.74,120 For the 6-membered ring in 12, a boat to chair conversion in moving to products released more strain than was alleviated in the 5-membered case. Magnus et al. found that the strain increased by 1.5 kcal·mol−1 (MM2/PRDDO) in moving from R to TS for 13, but that strain was relieved by 6.0 kcal·mol−1 for the six-membered ring compound 12.74,120 This release is attributed to a boat-chair strain release mechanism available to the 6-membered ring, but not the 5-membered ring. The crystal structure of 12 showed the 6-membered ring to be in a boat conformation (12b), which is higher in energy by 0.3 kcal·mol−1 as calculated by MI. While crystallization led to the boat conformation, our calculations suggest that it is likely that the chair conformation will co-exist with the boat conformation. For this reason the strain-release mechanism must be reinvestigated for a chair reactant, and not just a chair transition state. In comparing strain release factors between R and TS (ΔTS-Tot) it is clear that the same trends hold when considering a ground-state chair conformation of 12. As seen in Table IV and Table VII, strain release in moving towards the TS is −0.7, 1.5 and 6.3 kcal·mol−1 for 12, 14 and 13, respectively, mirroring the trend in activation free energies (all levels). Given the close agreement in calculated Rcd distances, the strain release mechanism appears to dominate control of the cyclization of fused ring system enediynes. The strain in the reactant conformation (MI calculations) which is comparable to the MM2 strain energy of the reactant, and that reported historically for enediynes (MM2 values) also predicts the correct ordering where higher strain energies in the reactant are indicative of higher strain release values in moving along the reaction coordinate.
Compound 14 is the sp3 version of 13 and cyclized faster than 13 because of a more favorable release of strain in the TS. Similarly, the oxidized version of 12, namely the hydroxyl compound 12c, reacted so quickly that reaction parameters could not be determined at 273 K. Partial characterization of the cyclization of 12c using methodology MI:Sum provided aΔH‡ of 19.1 kcal·mol−1 at 310 K, one of the lowest enthalpic barriers calculated in this study. The corresponding sum-corrected enthalpic barrier for 12 is 24.2 kcal·mol−1, and 18.9 kcal·mol−1 for the natural product-mimicking 15. For this reason it is not surprising that the rate of reaction of 12c was too fast for measurement and lends support to the low barrier calculated. Comparison of 12 and 12c using strain relief and critical bond distance values from Table VII indicate that the lowering of the barrier height for 12c is due to both a shortening of Rcd by 0.16 Å and increased strain release in moving from R to TS of −4.8 kcal·mol−1, in comparison to only −3.5 for the carbonyl version.
Direct comparison of the free energy barriers for 12–14 was made using temperatures of 310, 397 and 360 K. The only reported ΔS‡ was for 12, whose −7.33 cal·mol−1·K−1 value was reasonably reproduce using MI (-5.4 cal·mol−1·K−1). From Table I, the measured free energy barriers are 26.3, 32.0 and 27.4 kcal·mol−1 and using MI:Sum are 25.9, 31.4 and 26.3 kcal·mol−1, respectively. Not only is MI:Sum providing a correct qualitative ordering of the reactivity of 12–14, and 12c, the mean unsigned error associated with the calculated free energy barriers is only 0.5 kcal·mol−1 with the calculated values lying within the experimentally determined error bars.
Compounds 28 and 29 share a similar scaffold, however the bridgehead keto/secondary alcohol functionality has been moved in relation to the enediyne. Similar to compound 12, the secondary fused ring is 6-membered. The calculated Rcd distances for 28 and 29 are equivalent to 12, however less strain is released in moving to the transition state, and slight strain build up is evident (1.1 and 1.4 kcal·mol−1). According to MI-sum, at 310 K, compounds 28 and 29 cyclize with free energy barriers of 24.7 and 25.5 kcal·mol−1, respectively. Experimental results for 28 are not available, but the measured ΔG‡ of 26.5 kcal·mol−1 for 29 is in close agreement with the MI-sum value of 25.5 kcal·mol−1 (ΔH‡ in near quantitative agreement). Nearly identical calculated reaction parameters for 28 and 29 indicate that the pendant methoxyphenyl group plays little mediating role in the cyclization reaction, at most raising the barrier by 1 kcal·mol−1. Maier and Greiner attributed a lower calculated barrier to cyclization (2.9 kcal·mol−1 using a reparameterized form of MM2) for 28 versus 29 to the electronic effect of the substituted phenyl.120,121 This lowered barrier was supported by the inability to isolate the reactant 28 at room temperature, albeit the differential calculated here was only 1 kcal·mol−1.
The calculated strain values for 28 and 29 as compared to 12–14 show these two compounds to have strain build up lying in between 12 and 14, and the reported free energy barriers support this strain attenuation. One would expect differential strain release of 28 and 29 to be minimal, as is reflected in their calculated values (Table VII).
Molecular modeling of smaller enediyne systems serves to elucidate control factors in the cyclization reaction. Identification of acceptable methodology in modeling these reactions must allow the extrapolation to larger system size to be computationally tractable. Rigorous treatment at the CCSD(T) level of systems much larger than 1 becomes insurmountable. With the goal of modeling larger natural product-like species, we included three molecules that represent expanded reactant systems.
Two of these structures, 21 and 22 characterized by Nicolaou et al., were synthesized as model systems of the dynemicin core.122 These two structures were pursued since the non-aromatized version, 30, cyclized spontaneously. Using MMX, the authors determined the Rcd distances to be 3.13, 3.23 and 3.23 Å for 30, 21 and 23, respectively. Experimental ΔG‡ for 21 and 22 were reported as 22.6 and 25.7 kcal·mol−1. Since the Rcd values for 21–22 were below the cut-off value for spontaneous closure, Nicolaou et al. identified these compounds as novel exceptions to the rule. They explain that the Rcd rule only applies to strain energy in the reactant, and that differential strain in moving through the TS can lead to exceptions. The increased barrier is attributed to the relative gain in resonance stabilization in 21 (25 kcal·mol−1; 61–36 kcal·mol−1) and 22 (23 kcal·mol−1; 84–61 kcal·mol−1) compared with 30 (36 kcal·mol−1).123 Since compounds 21–22 do not gain as much resonance stabilization, Nicolaou et al. rationalize the products’ free energy for 21–22 to be higher than 30 and by invoking the Hammond postulate the cyclization barriers are increased. The authors do not find strain energy to differentiate the three reactivities, but instead use the bond order of the ene-bond where the bond decreases in length from 1.45 to 1.42 to 1.33 Å for model systems of 22, 21 and 30. The trend in bond order is correlated with the increasing enthalpy of cyclization of the three molecules and attributed to increased strain in the transition states.
When comparing MI and MI:Sum derived free energy barriers at 303K, it is apparent that sum correction of the terms moves the calculated values further from experiment compared with the tight qualitative and quantitative agreement for MI when considering 21 and 22. Both structures are rather similar, such that MI:Sum indicates 21 will cyclize with a barrier 1 kcal·mol−1 lower than 22. Strain parameters and the critical bond distance Rcd are nearly equivalent for these two compounds.
When investigating 21, 22 and 30, Nicolaou and co-workers determined 30 to react spontaneously such that no reported thermodynamic parameters exist. Owing to the lack of characterization, one would expect that the calculated cyclization barrier for 30 would be lower and this is indeed the case as MI:Sum free energies of cyclization at 310K decrease from 21.7 to 21.6 to 19.4 kcal·mol−1 for 22, 21 and 30. Strain parameters and the critical bond distance do not afford explanations as to why 30 has a much reduced barrier (reduced by a minimum of 2.2 kcal·mol−1). Benzannulation of 1 to form 2 showed a decrease in the experimental barrier, however these two species have calculated free energy barrier heights that are within 0.8 kcal·mol−1 of each other. When structure 4 is benzannulated to form 27, the calculated barrier increases, a finding reflected in the experimentally determined values. Benzannulation may be the driver for 21 and 22 having a larger barrier then 30, but this would be contrary to the effect of benzannulation in going from 1 to 2 and from 4 to 27.
Our model systems MTC (31) and MTC2 (15) are relatively large molecular subsets of the natural product esperamicin. Early work done on this reduced system supported the lowering of the reaction cyclization barrier upon triggering.88 That study was performed using the various computational methods on the smaller MTC system, which is MTC2 without two key hydroxyl groups (one alpha to the keto, and one alpha to the sulfur groups). Calculations using MI found the Rcd distance in the untriggered (lacking SSSCH3 group) and triggered versions of MTC to be 3.55 and 3.21 Å. Structures 15 and 15b represent the triggered and untriggered species for MTC2. We find a similar effect on Rcd for the expanded system (compared to MTC) where triggering of 15 closes this gap from 3.49 to 3.18 Å. The shorter Rcd distance upon triggering was reflected in a lowering of the ΔH‡ of cyclization by 23.4 kcal·mol−1 at MI for MTC (triggered barrier of 20.1 kcal·mol−1). In comparing the enthalpic barrier of MTC to that of 3, the presence of the lone hydroxyl group and the fused ring were identified as of the source for the lowering of the barrier from 26.0 to 20.1 kcal·mol−1. Using the ONIOM methodology, BD(T)/cc-pVDZ on the 6 carbon core reactive fragment showed the barrier to cyclization (of this small fragment only) to be nearly identical for the triggered and un-triggered versions of MTC (20.8 and 20.4 kcal·mol−1), while the inclusion of the rest of the molecule using the lower level ONIOM layer indicated a lowering of the barrier from 40.0 to 17.7 kcal·mol−1 upon triggering. This relative difference matches closely the effect calculated at MI. Here, we calculate the change in the uncorrected electronic energy barrier between 15 and 15b to be 17.5 kcal·mol−1. Strain contributions captured by the lower layer seemed to play a more important role than Rcd in the lowering of the barriers. Though the strain parameters were not quantified in our earlier work for MTC, we calculate strain release for 15 triggered (−2.1 kcal·mol−1) and strain build up for 15b untriggered (12.0 kcal·mol−1). Calculation of the cyclization barrier of 3 using BD(T) was consistent with the MTC results indicating that the enthalpic barrier of the central core atoms was 20.0 kcal·mol−1. Since Rcd for this species was 3.39 Å, strain release again appeared more dominant than internuclear distance.
The model compound MTC lacked two additional proximal hydroxyl groups present in the natural product, and we sought to include this functionalization in the present study. The resultant structure, MTC2 (15) represents esperamicin A1, and lacks only the nucleic acid recognition elements. While no direct experimental measurement exists for the determination of the enthalpic barrier of 15, a close comparison can be made to the experimentally determined free energy of cyclization for calicheamicin (20). De Voss et al. determined ΔH‡ at 262 K to be 19.3 kcal·mol−1, which is in close agreement with our MI:sum value for 15 of 19.0 kcal·mol−1 at 298 K. Comparison of 15 with the structures of calicheamicin and esperamicin A1 depicted in Figure 1 indicate that the deletion of the carbamate arm and sugar substituents does not influence the energetics of the reaction to any extent. Incorporation of additional hydroxyl groups in moving from MTC to MTC2 had no effect on the enthalpic barrier at MI (Table IX: 20.1 versus 20.5 kcal·mol−1).
Considering the near quantitative agreement in several of the reaction parameters calculated in reference to tabulated experimentally determined values, and the close agreement of the remainder of the structures studied, MI:Sum is confirmed as an accurate and efficient method for treating Bergman cyclizations of an expanded set of reactants. Several qualitative trends among compounds sub-families are also well reproduced. Unfortunately, the accurate estimation of reaction enthalpy is not obtained with the same methodology. Across 29 enthalpic and free energy barriers for which error functions can be calculated, MI:Sum performs well, with a MUE of only 2.8 kcal·mol−1 when compounds 5–8 are left out of the error statistics. The error for MII:Sum increases to 4.1 kcal·mol−1 over the same data set. Reaction enthalpies, considering only 1 and 2, are not estimated well with MI:Sum as the MUE increases to 8.7 kcal·mol−1, however, MII:Sum has an MUE of only 1.9 kcal·mol−1 for this same parameter. The MUE diagnostic of the prediction of reaction enthalpies runs over only two structures since this parameter has been measured for only compounds 1 and 2. For this reason, we choose to identify MII:Sum as a cost effective way of estimating the reaction enthalpies, but not as accurate for calculating forward reaction barriers compared with MI:Sum. That being said, MII:Sum estimates of ΔHr are explicitly needed to back out the retro-cyclization barriers. For this reason, we have tabulated a composite ΔHretro‡ in Table X so that the reverse reaction barriers can be weighed against the second reaction barrier, which must be crossed for successful hydrogen atom abstraction. This composite value is the MI:Sum forward enthalpic barrier minus the MII:Sum reaction enthalpy, both of which independently represent the most accurate estimations for each quantity.
Retro barriers are important in determining whether enediyne species will likely be able to proceed to hydrogen atom abstraction, or return back towards the retro-Bergman cyclization product. An example would be compounds 4 and 27, where it has been experimentally determined that benzannulation alters the kinetics such that 27 has a lowered barrier to retro-cyclization.124 Cyclization was determined to be rate limiting for 4, but compound 27 had a rate-limiting step of hydrogen atom abstraction, which was attributed to both a lower retro-barrier and higher abstraction barrier. Values in Table IV and Table X support these experimental observations, since the S–T gap increases for 27, and the retro-barrier drops by a factor of 2. This is similar to the case for 1 and 2 discussed above.
The barrier to hydrogen atom abstraction varies depending on substrate, but the barriers range from 9–14 kcal·mol−1.43,67,85,125 With this range in mind, comparing retro-barriers in Table X, it is evident that enediynes 9, 21–22 and 27 are the only species that would have retro-Bergman barriers low enough to compete with hydrogen atom abstraction. All other barriers are large enough for hydrogen abstraction to dominate.
Questions concerning the standardization of experimental conditions in measuring cyclization barriers, in addition to variations in solvent and temperature, help to explain why it has been difficult to adopt a single explanation for what is controlling the cyclization rates. Two main hypotheses have driven the quest for new enediynes, one being the critical internuclear distance, and the other being strain build up in the reactant and/or strain release to the transition state/products. Qualitative trends have emerged within closely related analog sets, but discrepancies often hinder the applicability of a global model for reactivity control. Clearly, utilizing MI:Sum to calculate reaction enthalpic barriers has been shown to have high qualitative if not near quantitative agreement with an error of 2.8–3.7 kcal·mol−1 (depending upon error subset), a factor just outside the range of experimental errors found in Table I. Correlations between geometrical parameters and the experimentally derived thermodynamic parameters are not statistically significant, and consistently have nearly no predictive or explanatory power for this diverse and large data set (Table VIII). Conversion of all measured values to free energies instead of enthalpies did not provide better trend analysis (Supporting Information Table S1). What is needed is a consistent data set that has no variation in how the compounds are treated. Since this is exactly what MI:Sum ΔH‡ represents, it is possible to query the calculated values for trends which may or may not exist for a series of geometric parameters representing the reactants and transition states. In addition, modeling of the data allows for the inclusion of more than one parameter in the fitting of the barrier heights, which affords the ability of including multiple competing yet synergistic explanations. Table VIII provides the correlation coefficients between the calculated and experimental barrier heights, as well as R2 values. Direct correlation between the measured and calculated barriers is low across the entire set of structures (for graphical representation, see Supporting Information). The highest correlation between the experimental data set and any calculated parameter is an inverse correlation to the total strain energy of the cyclized product (-0.41). Far from satisfactory, this offers no insight into designing new molecules. Contrary to this statistical insignificance, several highly correlated terms can be identified to influence the calculated barrier heights. Most notably, a correlation of 0.86 to the length Rcd and a correlation of 0.92 to ΔTS, or the strain released in moving from reactant to transition state (see graphs, Supporting Information). The tabulated R2 values support the correlation effects. In probing the calculated parameters for trends, it was determined that compound 26 was an outlier such that this compound was not included in the statistics. Inclusion of 26 serves to lower the significance and to alter any predictive model that might fall out of the analysis, the later point driving our decision to remove it from the training set.
Predictive models to the measured experimental values do not show acceptable statistical significance (owing to variations in experimental conditions and measurement techniques), but a partial least squares fitting of the calculated enthalpies to Rcd and ΔTS provides the following dependence:
where Rcd is measured in Å, ΔTS is the B3LYP derived strain release in kcal·mol−1, and ΔH‡ is the barrier height in kcal·mol−1. This line is fit to the calculated enthalpic barriers at 310 K for consistency, but a similar equation is arrived at when using the calculated barriers at the temperature relevant for each experimental determination. This indicates that the temperature variation for reported experimental measurements has less contribution to the low statistical significance of the equation fit to measured values as opposed to variations in measurement techniques. In addition, the trend line fit to calculated free energy barriers takes the same form and magnitude as that to the enthalpic barriers. Equation 5 relates a combination of the strain release and the critical distance Rcd to the enthalpic energy barrier. Both coefficients take positive values indicating that both an increase in strain release and an increase in Rcd raise the barrier to cyclization. This predictive model has a correlation coefficient of 0.97 and an R2 of 0.93 to the MI:Sum ΔH‡ values. It is apparent from the synergy of these two terms that strain release and acetylenic proximity are not competing, but complimentary explanations. In the two subsets of compounds shown in Table VIII, 1) benzannulated species, or 2) ring closure compounds, it is clear that strain and internuclear separation hold differing weights. Strain release plays a more dominant role in the benzannulated species, while a combination of strain release and Rcd are important drivers for the ring closure compounds.
Due to the low correlation of calculated parameters to the experimental values, it is not surprising that regression fits to the experimental barriers were not statistically significant. In fact, the R2 of the fit drops from 0.93 to 0.23 when comparing the model built to calculated enthalpies, and experimental enthalpies, respectively. Recently, Capitani, et al., used a similar approach to explain the variation in reactivity of endiynes.72 Using 10 enediynes, they also find that the critical distance Rcd does not explain reactivity differences alone. They arrive at a predictive model to reactivity based off a combination of activation energy and the critical distance, where the R2 of their model is 0.87, slightly less significant than that presented here. Their model relies on using an average value for the critical distance at the transition state, while our model explicitly accounts for variation in strain energy for each individual transition state.
Among some of the additional questions surrounding the present work is the effect of solvent on the reaction barriers. It is apparent from the experimental data that large variation is possible in the measured thermodynamic parameters across a series of solvents including the question as to whether aromatic solvents interact in any specific manner. We think a better understanding of the factors governing cyclization rates has been established through the use of gas-phase calculations. Investigation of implicit and/or explicit solvation is beyond the scope of this study, though we expect that it should yield interesting insights into this class of compounds. Additionally, several questions remain to be answered about competing reactions given varying concentrations of reaction co-factors such as 1,4-cyclohexadiene. While we have invoked the effects low concentrations of 1,4-CHD have on the reaction kinetics, it is of future interest to more fully characterize competing reactions that may help to rationalize discrepancies in reported thermodynamic parameters discussed above.
We thank Prof. Chris Cramer (U. Minnesota) for his insight into the sum method correction, and reviewers for useful insight and critical reading of the manuscript. Acknowledgment is made to Research Corporation, to NIH, and to NSF for support of this work. This project was supported in part by NIH grant 1R15CA115524-01, NSF grant CHE-0457275, and by NSF grants CHE-0116435 and CHE-0521063 as part of the MERCURY high-performance computer consortium (http://mercury.chem.hamilton.edu).
Supporting Information Available: Additional information includes: Panel 1: Critical bond distance vs. calculated enthalpy; Panel 2: Calculated vs. experimental enthalpy at appropriate temp; Panel 3: Strain release from R to TS calculated with DFT vs. calculated enthalpy; Panel 4: Predicted vs. calculated enthalpy based on two parameter model; and Table S1: Enthalpies and Free Energies derived from Calculated Entropy (kcals·mol−1). Cartesian coordinates and electronic energies for stationary points are provided. This information is available free of charge via the Internet at http://pubs.acs.org.