How Robust Is the Reversible Steric Shielding Strategy for Photoswitchable Organocatalysts?

A highly appealing strategy to modulate a catalyst's activity and/or selectivity in a dynamic and noninvasive way is to incorporate a photoresponsive unit into a catalytically competent molecule. However, the description of the photoinduced conformational or structural changes that alter the catalyst's intrinsic reactivity is often reduced to a handful of intuitive static representations, which can struggle to capture the complexity of flexible organocatalysts. Here, we show how a comprehensive exploration of the free energy landscape of N-alkylated azobenzene-tethered piperidine catalysts is essential to unravel the conformational characteristics of each configurational state and explain the experimentally observed reactivity trends. Mapping the catalysts' conformational space highlights the existence of false ON or OFF states that lower their switching ability. Our findings expose the challenges associated with the realization of a reversible steric shielding for the photocontrol of Brønsted basicity of piperidine photoswitchable organocatalysts.

Introduction

Inspired by the regulation of the activity of enzymes through trigger-induced effects,[1] chemists are seeking to develop synthetic catalysts whose reactivity or principle of asymmetric induction can be modulated by external stimuli, such as temperature, pressure, or pH.[2] Light is another particularly attractive stimulus due to its noninvasiveness, high spatial and temporal resolution, and possibility of being manipulated precisely with modern optics.[3] In photoswitchable organocatalysis, electromagnetic irradiation induces a reversible transformation of a photochromic moiety incorporated into an organocatalytic system (e.g., (E) → (Z) isomerization of a double bond).[4] The photochrome must be carefully chosen to ensure that the photoinduced reaction occurs with a high efficiency in both the forward and reverse directions, meaning high concentrations of either configuration can be obtained. This transformation must then cause a significant alteration of the steric or electronic properties of the catalyst, which leads to the two states having intrinsically different reactivities.[5] The configurational state displaying a higher reactivity is termed the ON state, while the one with a lower reactivity is called the OFF state.

Different strategies have been developed to attain this modulation of chemical reactivity with light, including activation,[6] inhibition,[7] templation,[8−19] electronic alteration,[20−28] and reversible steric shielding. The latter is a classic, conceptually elegant and, in principle, straightforward design strategy, exploited as early as 1981,[29,30] whereby access to either the substrate-binding site, or the catalytically active site, is restricted by a blocking group in the OFF state (Scheme ). In the ON state, the photoinduced transformation/isomerization allows the blocking group to be moved away from the active site and the molecule to become catalytically active.[31] In an ideal scenario, there should be a one-to-one correspondence between a catalyst’s configurational state (e.g., (E)- or (Z)-isomer) and it being in either the ON or OFF state (i.e., active-site shielded or deshielded), meaning that whenever the catalyst is, for example, in the (Z)-configuration, it is also in the ON state (while the (E)-configuration is always OFF). Unfortunately, the conformational behavior of a photoswitchable organocatalyst can complicate this picture and break this correspondence down: since organocatalysts are generally rather flexible molecules with many low-energy, thermally accessible conformations, the same configuration might display both ON and OFF activity.

Scheme 1

Photoswitchable Azobenzene-Based Piperidine Organocatalysts for the Henry Reaction Exploiting Reversible Steric Shielding

In the (E)-configuration, a steric shield blocks access to the active site. In the (Z)-configuration, photoinduced switching frees access to the active site. This strategy relies on a one-to-one correspondence between the configuration and the ON/OFF state.

To understand the catalytic behavior of flexible organocatalysts that have been rationally designed based on reversible steric shielding (but also on hydrogen bonding templates, π-stacking, anion-π interactions, etc.),[32−36] the vast conformational space associated with their rich morphology must be thoroughly explored. Although conformational analysis is a familiar topic to organic chemists,[37] computational studies on organic systems are frequently limited to static density functional theory (DFT) depictions of isolated minima, aligning with textbook two-dimensional presentations of reaction mechanisms. This approach can fail, sometimes even leading to erroneous conclusions,[38] when medium-sized or large, highly flexible molecules that evolve on complex free energy surfaces (FES) are investigated.[39−43] While the conformational sampling of medium-sized organic systems can be performed with tools such as CREST,[44,45] the latter aims at improving the description of ensemble averages through Boltzmann reweighting of the optimized structures rather than at providing access to the free energy landscape. Ab initio molecular dynamics is limited to short time intervals (∼102 to 103 ps) in which chemically relevant interconversions between important states might not be explored.[46−49] Using potentials derived from semiempirical methods, such as density functional tight binding (DFTB)[50−56] or machine-learned potentials[57−63] that can reach the accuracy of quantum mechanical computations at a fraction of their computational cost, allows the simulation times to be lengthened, yet owing to the complexity of the free energy landscape of photoswitchable organocatalysts, visiting all possible conformational/configurational regions might still be impossible. Recently, our group exploited enhanced sampling techniques, combined with low-cost electronic structure computations[64] or neural network[65] and kernel-based potentials,[66] to address organic chemistry problems connected to fluxional molecules, such as organocatalysts[67,68] and molecular rotors.[69] Within this context, Hamiltonian or temperature replica exchange molecular dynamics (REMD)[70] are well-suited sampling approaches because they do not require any a priori knowledge of the dominant conformational regions or of the relevant collective variables characterizing them and are able to generate unbiased canonical sampling at a particular temperature.[66]

In this work, we use parallel tempering REMD simulations at the DFTB3 level[71−73] to investigate the FES and the corresponding configurational and conformational behavior of an early example of photoreversible steric shielding, namely the azobenzene-tethered N-alkylated piperidine base catalysts developed by Hecht and co-workers (Scheme ).[74−77] This elegantly designed system consists of a piperidine Brønsted base and a photochromic azobenzene moiety that are rigidly and orthogonally positioned through the spiro junction of an isobenzofuranone ring. In the thermodynamically more stable (E)-isomer, the 3,5-substituents on the phenylazo unit shields the basic site; UV irradiation affords the metastable (Z)-isomer and exposes the nitrogen lone pair, enhancing the basicity by almost one order of magnitude and the rate of the Henry reaction[78,79] between nitroethane and p-nitrobenzaldehyde. Analysis of the potential energy surface of PSa–c, obtained through the computation of a handful of static DFT geometries, in combination with NMR experiments, was previously used to optimize the catalyst’s structure.[75] However, this approach cannot fully explain the observed reactivity trends of photoswitches PSa–c (Scheme ). Our simulations shed light on the subtle interplay between conformational states and activity and expose the possible drawbacks of reversible steric shielding as a design strategy to achieve light-gated modulation of chemical reactivity.

Computational Details

All potential energy computations were performed in the gas phase at the DFTB3/3OB[72,73] level in combination with the D3BJ[80] dispersion correction, as implemented in DFTB+ software.[81] DFTB is two to three orders of magnitude faster than ab initio and DFT methods, making it particularly attractive in applications to large molecules and condensed phase systems, especially when extensive sampling is important to the reactive process of interest. Static DFT and DFTB3 computations were performed with Gaussian16[82] and AMS2020.1[83,84] software packages to assess the impact of the level of theory on the barrier heights and relative stability of different conformations. Parallel tempering (PT) replica exchange simulations were performed using the REMD@DFTB3[64] protocol implemented in i-PI[85] to sample the canonical (i.e., NVT) ensemble at the target temperature of 300 K for each photoswitchable organocatalyst. In PT-REMD, a series of energetically independent trajectories (called replicas) of the same system at different temperatures is simulated, allowing complete configurations from replicas at these temperatures to be occasionally exchanged. The exchange between replicas propagated at higher temperatures, where energy barriers can be overcome and a large amount of the FES is explored, and those at a lower temperature, constrained to local minima of the FES, ensures converged statistical sampling at a target temperature and thorough exploration of the free energy landscape. The (E)- and (Z)-configurations of each species were sampled in independent simulations, using a harmonic restraint to maintain the organocatalyst in a specific configuration throughout the simulation. The restraint was introduced in the dynamics using PLUMED[86−88] software. The simulations included 16 replicas with temperatures ranging from 300 to 1200 K separated by logarithmic intervals. A time step of 0.6 fs was used in the dynamics, with a Langevin thermostat to control the temperature. To ensure the convergence of the sampling, we ran the simulations for 1.2 ns and analyzed the evolution of the integrated free energy differences between basins at 300 K. As the crossing of energy barriers is associated with the slowest dynamic modes of the systems, their convergence is a reliable estimation of the overall convergence of the sampling (see the Supporting Information). The interactive plots were constructed with Python framework Dash for web applications (https://plotly.com/dash/).

While it is well-known that dispersion-corrected DFTB3/3OB leads to potentials significantly flatter[64,66,72,89,90] than those computed with a higher electronic structure level, the consequences of these shortcomings have been extensively discussed.[64,66,72,90] The reader is thus pointed to refs (64), (66), (72), (89), and (90) for further technical considerations on the topic. Importantly, it has been demonstrated by us and others that the low-cost computational level provides a reliable description of organic chemical reactions[90] and is capable of identifying all the relevant conformational regions accessible around 300 K. For the photoswitchable systems investigated herein, differences between static scans of selected collective variables at the DFTB3 and DFT levels (i.e., without conformational entropy corrections) and comparisons with the FESs (accounting for conformational entropy) are discussed in the Supporting Information. In addition, the consequence of the relative flatness of the low-cost potential on the relative population of catalytically active or inactive states by PSa–c is discussed. Overall, it is confirmed that accounting for the full conformational entropy and anharmonic effects is more transformative for the relative population of the conformational regions than the use of higher-level potentials.

Results and Discussion

Analysis of the Free-Energy Maps

The conformational analysis of PSa–c has previously[75] been reduced to the limited number of chemically intuitive static representations 1–5 shown in Scheme , and three pathways have been considered to account for the catalysts’ residual activity in the (E)-state: (1) inversion of the piperidine N-atom; (2) inversion of the piperidine chair, that is, “ring flip”, leading to the more stable conformer 4 after N-inversion; (3) rotation of the blocking group around the (N=N)–Ph bond, placing the benzofuranone and phenylazo groups in an orthogonal orientation and allowing the substrate to be deprotonated by piperidine in its most stable chair conformation with the R1 group in the equatorial position. Similar processes (N-inversion, ring flip, and rotation) can also occur in the (Z)-configuration but are expected to be less impactful on the catalytic activity as the N-lone pair should always remain accessible. On the other hand, a finite temperature exploration of a molecule’s free energy landscape (Figure ) may alter, if not even reverse, pictures provided by static relative energy computations.[64]

Scheme 2

Schematic Representation of Local Minima and Transition States Rationalized according to Chemical Intuition and Static DFT Computations[75]

Reproduced from Stoll, R. S.; Peters, M. V.; Kuhn, A.; Heiles, S.; Goddard, R.; Bühl, M.; Thiele, C. M.; Hecht, S. J. Am. Chem. Soc. 2009, 131 (1), 357–367. Copyright 2009 American Chemical Society.

Figure 1

Two-dimensional representation of the free energy landscape corresponding to the NVT ensemble at 300 K obtained from the REMD@DFTB3 of (E)- (top) and (Z)-PSa–c (bottom). The relevant collective variables are shown at the top. The regions of the landscapes are colored in green or red according to the ON or OFF function of the photoswitches. The orange and blue points represent the local minima and transition states (labeled according to Scheme ) optimized at the DFTB3 level from the static geometries reported in ref (75). Distribution histograms according to d and θ are shown on the top and right axes of the plots, respectively.

Figure shows the 2D FES corresponding to the NVT ensemble at 300 K of PSa–c represented in terms of two collective variables, d and θ. The former is the shortest distance between the nitrogen atom of the piperidine ring and any other atom of the phenylazo group. The latter is the distance between the piperidine’s N-atom and the plane passing through the three carbon atoms the nitrogen is bound to (R1-C and two C-atoms of the piperidine ring). If the nitrogen is located below this plane, θ adopts negative values; if N is above the plane, θ is positive (see Figure ). The conformational region with θ < 0 is labeled I, and the one with θ > 0 is labeled II. Each FES is divided into two areas, colored green or red, according to the values of d and θ. The red area indicates that the catalyst’s basic site (i.e., the nitrogen’s lone pair) is sterically shielded, while the green one indicates that it is deshielded.[91] Interactive plots showing the conformational states used to construct the FES are freely available as a Heroku deployed app (https://photoswitch-exp.herokuapp.com/).

In the (E)-state plots, region I is the conformational region associated with 1 and 3 and the transition states for ring flip and rotation; region II contains 2 and 4, while transitions between the two regions are associated with the inversion transition state (Scheme ). For comparison with the static picture, the structures previously determined with DFT computations,[75] reoptimized at the DFTB3 level, are included as the orange and blue points in Figure .

Region I of (E)-PSa has the largest conformational entropic contributions, as indicated by the larger basin and small REMD free energy differences between the ring-flipped conformers (the DFTB3 free energy difference between the red and green areas of basin I of (E)-PSa is 1.3 kcal/mol, vs 2.1 kcal/mol for (E)-PSc). Since there are many possible twist or twist-boat intermediates and transition states, depending on which C–C or C–N bond of the piperidine ring is rotated relative to its opposite counterpart and in what direction, the generation of ring-flipped conformers cannot be thought of as an individual process but rather as the population of a continuum of twist-boat (or twist) conformational states, which cannot easily be captured by a static minimum-energy structure search approach. Compared to (E)-PSb–c, the broader conformational region I of (E)-PSa is associated with the population of false OFF states in which the benzofuranone and phenylazo rings are oriented orthogonally to each other through rotation around the (N=N)–Ph bond (see, e.g., 6 on the bottom left corner of Figure ). In (E)-PSb and PSc, the bulkier R2 groups (Bu or 2,6-Me2Ph) help more effectively shield the lone pair, making region I more structured in the red area of the plot.

The small energy difference between region I and II of (E)-PSa (0.5 kcal/mol) suggests that N-inverted species with the R1 (Me) group in the axial position and rearranged conformers with equatorial R1 and exposed N-lone pairs are significantly populated. Experimentally, the 13C NMR shift of the N-methyl carbon of (E)-PSa was found to lie in between the shifts for axial and equatorial positions.[75] This suggests that N-inversion is fast on the NMR timescale, in agreement with the significant population of region II and fast exchange of the N-methyl between axial and equatorial positions. In PSb and PSc, basin II is smaller and possesses higher relative free energies (1.2 and 1.7 kcal/mol, respectively), indicating that N-inversion in these catalysts is suppressed. This relates well with the larger size of the R1 substituent (Bu group) and its unlikeliness of adopting the axial position in a six-membered ring.

The conformational space in the (Z)-configuration is narrower (the FES possesses a more distinct minimum) than in the (E)-state (the FES is flatter). The minimum of region I of (Z)-PSa is clearly located within the green area of the plot. Owing to the small ratio between the van der Waals volumes of R1 (Me, 25.8 Å3) and R2 (Bu, 72 Å3),[92] only a small fraction of false ON states in which the active site is shielded is populated. Like in the (E)-state, regions I and II of (Z)-PSa possess similar free energies and an inversion at the piperidine’s N-atom with an exchange of the axial lone pair and equatorial methyl group has a high probability of occurring. Region II of (Z)-PSb–c is less frequently visited due to the unlikeness of inverting the bulkier R1 (Bu) group, which would lead to significantly unfavorable 1,3-diaxial interactions. On the other hand, (Z)-PSb and (Z)-PSc have a higher probability than (Z)-PSa of visiting the red area of the plot; despite having the same R2 group as PSa (Bu), PSb is more likely to populate false ON states in which the blocking group shields the catalytically active site due to a rotation of the benzofuranone–azo bond (e.g., structure 7 on the bottom right corner of Figure ). Although the meta-xylene R2 group of PSc (110.6 Å3) is bigger than tert-butyl in PSb (72 Å3),[93] the aromatic ring is “flatter” than the alkyl substituent and, therefore, less likely to shield the basic site in the (Z)-configuration. Thus, PSc has a lower probability than PSb to populate false ON states.

Overall, an inspection of the Figure plots shows that none of the catalysts examined is ideally tuned for reversible steric shielding: in an ideal photoswitchable catalyst, the conformational space of each configuration should be totally restricted to either the red or green area. However, within each (E)- or (Z)-plot, basins located both in the green and red areas are populated. The broad conformational space of (E)-PSa and (Z)-PSb–c is illustrative of the importance of capturing the full anharmonic and entropic contributions to the free energy landscape: computing a “static” population of few representative conformers cannot properly account for the existence of false OFF and false ON states that reduce the catalysts’ switching ability.

Relating Conformational Behavior with Experimental Reactivity

Kinetic studies on the catalytic performance of the three photoswitchable catalysts showed that, in the (E)-configuration, PSa is the most active, while PSc leads to an essentially trace yield of β-nitro alcohol (Table ).[75] This is consistent with the free energy maps in Figure , which show that the Bu group of PSc–b serves as an efficient conformational anchor for the six-membered piperidine ring and that the meta-xylene R2 group improves the efficiency of blocking the active site, localizing the minimum of the FES of (E)-PSc in the red inactive region. On the other hand, (E)-PSa can access enthalpically and entropically favored conformational states in the green active region due to inversion at the piperidine’s N-atom and rotation around the (N=N)–Ph bond.

Table 1

Photochemical, Kinetic, and Thermodynamic Experimental Data of Piperidine Bases PSa–c Taken from Ref (75)a

	PSa	PSb	PSc
PSS (Z/E)b	90:10	90:10	>90:10
t1/2c[h]	268	286	466
kEd [10–6 s–1]	5.0	0.96	0.39
YEe (%)	25	6	0
kZf [10–6 s–1]	22.0	13.0	14.0
YZg (%)	72	47	57
krel (kZ/kE)	4.4	13.5	35.9
ΔpKah		0.8	0.7

Reproduced from Stoll, R. S.; Peters, M. V.; Kuhn, A.; Heiles, S.; Goddard, R.; Bühl, M.; Thiele, C. M.; Hecht, S. J. Am. Chem. Soc. 2009, 131 (1), 357–367. Copyright 2009 American Chemical Society.

Photostationary state (PSS) obtained by irradiation at 365 nm.

Half-life of the (Z)-isomer at 20 °C.

Rate constant of the Henry reaction using puI(E)-isomer.

2-Nitro-1-(4-nitrophenyl)propan-1-ol yield usingIre (E)-isomer.

Rate constant of the Henry reaction extrapolated to 100% (Z)-isomer.

2-Nitro-1-(4-nitrophenyl)propan-1-ol yield extrapolated to 100% (Z)-isomer.

Difference of pKa values, i.e., pKa(PSS) – pKa(E), obtained from titration with triflic acid using neutral red as a reference base.

The free energy plots of the (Z)-configuration suggest that PSa has a lower probability than PSb or PSc of visiting the red area of the FES thanks to the small size ratio of its R1 and R2 groups and their inability to effectively shield the lone pair. Consistently, (Z)-PSa displays the highest rate and yield in the Henry reaction (Table ). Owing to region I of their FES being broader, along with the higher energy difference between basin I and II, (Z)-PSb and (Z)-PSc can access inactive states through rotation around the isobenzofuranone core and thus have a lower catalytic activity.

The switching ability of catalysts PSa–c can also be rationalized according to the conceptual plot shown in Figure , where the relative rate (k/k) is expressed as a function of the conformational behavior. At the top of this plot lies the ideal photoswitchable organocatalyst, whose distribution of conformations in the (E)-state is always inactive, while in the (Z)-state is always active. This would correspond to an ideal free energy map in which the conformational space of each configuration is totally restricted to either the red or green area. At the bottom of the left and right sides lie the “worst” catalysts, namely those in which false OFF states are significantly populated (left side) or true ON states are scarcely populated (right side). On traversing the plot left to right, populating false OFF states becomes increasingly more difficult due to a more structured (E)-state and k decreases, while the population of true ON states (i.e., k) is assumed to remain approximately constant until the peak is reached. Further progression to the right corresponds to a reduction of k (with k ca. the same) since the (Z)-configuration becomes more flexible, leading to false ON states being populated.

Figure 2

Conceptual plot showing the correlation between the conformational behavior of photoswitches PSa–c and their switching ability (kreli.e., k/k). RSS is a hypothetical reversible steric shielding descriptor.

PSa lies toward the bottom of the left side of the plot due to its highly flexible (E)-configuration and significant population of N-inverted and (N=N)–Ph rotated conformers; PSb is located on the opposite side, because of the likeliness of populating (Z)-conformational states in which the phenylazo unit shields the piperidine’s N-atom. Interestingly, PSc, although possessing the highest krel (Table ), is not located at the top of the plot because, compared to PSa, it has a less structured (Z)-state, and true ON states are populated to a lesser extent, implying the possibility of further improvement. Reversible steric shielding is clearly a double-edged sword: while the introduction of the large tert-butyl group on the piperidine ring and of the meta-xylene substituent on the phenylazo unit effectively maximizes the population of true OFF states in the (E)-configuration, locking the piperidine ring with an axial lone pair adequately shielded by the R2 substituent, it simultaneously leads to a less structured (Z)-configuration in which false ON states are more frequently visited.

Although the present analysis was performed in the gas phase (the inclusion of an explicit solvent is not easily compatible with PT)[94] and at a fairly low electronic structure level (REMD@DFTB3), the key findings related to the conformational behavior of the photoswitchable catalysts examined remain valid. In particular, solvent effects tend to flatten the DFT (or higher)-based FES (see Table S2).[66] Overall, the relative population of catalytically active or inactive states of PSa–c, corroborated by comparison to the experimentally observed reactivity trends, is not expected to significantly change. While obtaining both highly accurate energetics and converged statistical sampling has typically hampered the appropriate description of the flexibility of organic molecules, efforts in our group are being made to achieve ab initio accuracy by correcting semiempirical potentials with machine learning models, including implicit[66] or explicit treatment of solvation.[65,95]

Conclusions

In this work, we have illustrated the importance of thoroughly exploring of the free energy landscape of flexible photoswitchable organocatalysts to understand their catalytic performance. Without any a priori knowledge of the system, the complex FES of three N-alkylated azobenzene-tethered piperidine photoswitches were successfully mapped, the energetic basins corresponding to the structural minima identified, and the experimentally observed reactivity trends rationalized according to the catalysts’ conformational behavior. Our findings show that the lack of highly structured configurations and the access to entropically favored regions in one state can limit the effectiveness of reversible steric shielding as a design strategy to attain photoregulation of chemical reactivity. While large substituents and steric clashes may have a beneficial effect in one configuration, leading to a limited number of false OFF states being populated, they can simultaneously result in the population of false ON states in the other, eroding a catalyst’s photoswitching ability. This example explicitly illustrates the difficulty of achieving structure-based optimization on the basis of a static picture. In this respect, free energy sampling appears necessary to predict the impact of the floppy or bulky nature of different substituents on the conformational landscape of the ON/OFF states in the search of improved azobenzene derivatives with a higher structural rigidity.